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Self-consistent  models  of  barred  spiral  galaxies  based  on  the  observed  properties 
of  NGC  3992,  NGC  1073,  and  NGC  1398  are  constructed  and  analyzed.  The  method  of 
model  construction  is  a slight  modification  of  the  technique  developed  by  Contopoulos 
and  Grosbpl  for  the  case  of  unbarred  spirals.  The  main  factors  which  influence  self- 
consistency  are  the  amplitude,  pitch  angle,  scale  length  and  z-thickness  of  the  spirals,  the 
mass  of  the  bar,  the  angular  velocity  of  the  bar/spiral  pattern,  the  central  surface  density 
and  scale  length  of  the  disk,  and  the  central  value  and  slope  of  the  velocity  dispersion. 

Stochastic  orbits  whose  Jacobi  constants  lie  between  the  values  at  the  Lagrange 
points  L\  and  L4  are  found  to  play  a significant  role  in  supporting  self-consistent  spiral 
structure,  especially  in  the  regions  just  beyond  the  ends  of  the  bar.  Stochastic  orbits 
whose  Jacobi  constants  lie  below  this  interval  tend  to  fill  more  or  less  uniformly  either 
rings  in  the  outer  disk  or  ovals  in  the  bar  region,  depending  on  the  regions  to  which  they 
are  confined.  Stochastic  orbits  whose  Jacobi  constants  lie  above  that  of  L4  also  tend 
not  to  support  any  imposed  structure.  The  model  bars  are  predominantly  comprised  of 
elongated  orbits  trapped  around  the  X\  family  and  terminate  close  to  corotation. 
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CHAPTER  1 
INTRODUCTION 


Spiral  structure  in  galaxies  was  first  observed  by  William  Parsons,  Third  Earl  of 
Rosse,  in  several  nebulae  that  had  been  catalogued  earlier  by  William  Herschel.  Because 
this  discovery,  made  in  the  mid-nineteenth  century,  preceded  the  establishment  of  these 
nebulae  as  separate  galaxies  external  to  the  Milky  Way,  they  became  known  simply  as 
spiral  nebulae.  In  addition,  their  spiral  shapes  seemed  to  suggest  that  the  systems  rotate 
about  an  axis  perpendicular  to  the  plane  of  the  spiral.  Confirmation  of  this  hypothesis 
was  provided  in  1914  by  Vesto  Slipher,  whose  spectroscopic  observations  of  a number 
of  these  spiral  nebulae  revealed  the  Doppler  shifts  produced  by  rotation.  The  nature  of 
the  spiral  nebulae,  one  of  the  subjects  of  the  so-called  “great  debate”  between  Huber 
D.  Curtis  and  Harlow  Shapley  at  the  National  Academy  of  Sciences  in  April  1920,  was 
firmly  established  in  the  early  1920s  by  Edwin  Hubble. 

In  1923,  Hubble,  using  the  100-inch  telescope  at  Mount  Wilson,  was  able  to 
resolve  the  outer  parts  of  the  spiral  nebulae  M3 1 and  M33  into  multitudes  of  apparent 
point  sources.  Although  at  first  unable  to  determine  whether  these  point  sources  were 
individual  stars,  he  soon  found  that  some  were  indeed  Cepheid  variables.  Using  the 
Cepheid  period-luminosity  relation,  originally  discovered  in  1912  by  Henrietta  Leavitt, 
Hubble  was  able  to  show  conclusively  that  M31  and  M33  are  at  very  large  distances, 
and  therefore  galaxies  external  to  the  Milky  Way. 
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In  1926,  Bertil  Lindblad,  working  within  Shapley’s  concept  of  the  Milky  Way 
(i.e.  that  the  Sun’s  position  is  generally  in  the  galactic  plane  but  relatively  far  from 
the  galactic  center),  developed  a mathematical  model  for  its  rotation.  This  model 
depicted  the  galaxy  as  consisting  of  a number  of  subsystems,  each  rotating  with  its  own 
characteristic  angular  velocity  and  degree  of  flattening  about  a common  axis.  Jan  Oort, 
in  1927  and  1928,  demonstrated  the  basic  correctness  of  Lindblad’s  model  and  showed 
that  the  galactic  disk  is  in  a state  of  dijferential  rotation.  That  is,  stars  closer  to  the 
galactic  center  orbit  around  it  with  higher  angular  velocities  than  those  farther  out  in 
the  disk.  At  about  the  same  time  Lindblad  began  to  try  to  understand  the  nature  of 
spiral  structure,  a problem  he  worked  on  until  his  death  in  1965. 

Perhaps  the  simplest  interpretation  of  galactic  spiral  structure  is  that  it  represents  a 
pattern  of  material  arms  that  rotate  as  a rigid  body  about  the  galactic  center.  Differential 
rotation,  however,  presents  this  interpretation  with  a serious  difficulty  called  the  winding 
dilemma.  Consider,  for  example,  a galactic  spiral  arm  of  uniform  pitch  angle  I'o  = -30° 
at  time  r = 0.  Also,  suppose  that  the  galaxy’s  angular  rotation  rate  fl  depends  on  radius 
r (i.e.  the  galactic  disk  rotates  differentially).  The  azimuth  of  the  arm  at  radius  r and 
time  t is  given  by 


(1-1) 


The  pitch  angle  i at  radius  r and  time  t is  defined  by 


(1-2) 


For  a typical  galaxy  with  a flat  rotation  curve,  Q,{r)r  ^ 200  km  s \ r w 10  kpc,  and 
t PS  10^®  yr.  In  this  time  the  pitch  angle  i would  have  decreased  from  the  original  -30° 
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to  a value  of  -0.°28.  Therefore,  any  material  arms  would  have  wound  up  and  rendered 
the  spiral  structure  unrecognizable  on  timescales  comparable  to  the  age  of  the  galaxy. 

Binney  and  Tremaine  (1987)  summarize  four  possible  solutions  to  the  winding 
dilemma.  Firstly,  spiral  structure  may  represent  a statistical  equilibrium  in  which 
the  age  of  any  given  arm  is  quite  young.  The  idea  here  is  that  clumpy  features  are 
continuously  produced  in  the  galactic  disk  and  are  sheared  off  into  spiral  segments 
by  differential  rotation.  A spiral  arm,  then,  would  simply  represent  an  aggregation  of 
these  shorter  spiral  segments.  While  this  chaotic  spiral  arm  hypothesis  is  believed  to 
be  applicable  to  flocculent  galaxies  (e.g.  Toomre  1989),  it  has  difficulty  accounting  for 
the  striking  global  coherence  exhibited  by  the  spiral  patterns  of  grand  design  galaxies. 
Secondly,  spiral  structure  may  be  the  temporary  result  of  a recent  tidal  interaction  with 
a companion  galaxy.  In  this  case  the  companion  excites  a transient  global  spiral  density 
wave  in  the  primary  galaxy.  The  nature  of  density  waves  is  discussed  in  more  detail 
below.  Models  of  this  type  are  able  to  reproduce  many  features  observed  in  grand  design 
galaxies,  such  as  dust  lane  location  and  strength  as  well  as  the  main  properties  of  the 
neutral  hydrogen  distribution  and  velocity  field.  According  to  Binney  and  Tremaine 
(1987,  p.  394),  however,  “they  cannot  account  for  all  spirals  because  encounters  with 
massive  companion  galaxies  in  favorable  orbits  are  not  common  enough.”  Thirdly,  spiral 
structure  may  represent  a detonation  wave  of  star  formation,  driven  by  supernovae 
explosions  or  expanding  HII  regions,  that  propagates  around  the  disk.  The  wave  is 
sheared  into  a trailing  spiral  by  differential  rotation  and  ultimately  settles  down  with  a 
fixed  shape  and  pattern  speed.  This  self-propagating  star  formation  hypothesis,  though, 
requires  finely  tuned  star  formation  rates  and  does  not  adequately  address  the  problems 
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of  the  broad  arms  of  the  older  disk  components,  the  very  regular  patterns  seen  in  some 
grand  design  galaxies,  and  the  streaming  observed  along  spiral  arms.  Finally,  spiral 
structure  may  represent  underlying  waves  in  the  density  and  gravitational  potential  of  the 
disk  itself.  The  orbiting  stars  and  gas  adjust  their  motions  such  that  they  tend  to  linger 
near  the  minima  of  the  potential;  moreover,  the  gas  in  these  regions  is  compressed, 
thereby  inducing  rapid  star  formation  that  generates  the  bright  young  stars  and  HII 
regions  which  delineate  the  visible  spiral  pattern.  While  the  material  that  constitutes 
the  density  waves  is  constantly  changing,  the  waves  themselves  are  neutrally  stable 
and  represent  a quasi-stationary  spiral  pattern  that  rotates  as  a rigid  body  with  angular 
velocity  Q.p.  Bertil  Lindblad  (1961,  1963)  originally  proposed  such  a scenario,  but  his 
approach,  which  emphasized  the  role  of  individual  stellar  orbits,  was  not  well-suited 
for  quantitative  global  analysis.  C.  C.  Lin  and  Frank  Shu  (1964)  first  provided  the 
mathematical  formalism  of  the  spiral  density-wave  theory.  Therefore,  the  idea  that 
spiral  structure  represents  a quasi-stationary  density-wave  has  become  known  as  the 
Lin-Shu  hypothesis.  Since  the  present  thesis  concerns  the  structure  of  barred  spiral 
galaxies,  we  must  also  consider  those  phenomena  associated  with  the  presence  of  a bar. 

Bars  are  found  in  a large  fraction  of  disk  galaxies.  In  fact,  depending  upon  the 
reference  catalog  cited,  between  25%  and  35%  of  all  disk  galaxies  are  classified  as 
strongly  barred  (Sandage  and  Tammann  1981;  de  Vaucouleurs  et  al.  1976;  Nilson 
1973).  An  additional  fraction,  about  whose  size  there  is  greater  disagreement,  is 
classified  as  possessing  weak  bars  or  oval  distortions.  The  class  of  barred  galaxies 
is  very  heterogeneous.  In  it  we  find  galaxies  that  span  the  entire  range  of  Hubble  types. 
In  addition,  “other  properties,  such  as  the  size  of  the  bar  relative  to  the  host  galaxy,  the 
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degree  of  overall  symmetry,  the  existence  of  rings  and  numbers  (and  position  relative  to 
the  bar)  of  spiral  arms  in  the  outer  disc,  the  gas  and  dust  content,  etc.,  vary  considerably 
from  galaxy  to  galaxy”  (Sellwood  and  Wilkinson  1992,  p.  1). 

Observations  of  bars  generally  concern  either  the  light  distribution  or  the  velocity 
field.  Analyses  of  the  light  distributions  reveal  that  bar  radii  (i.e.  semimajor  axes)  are 
always  less  than  the  disk  radii  (defined,  for  example,  by  the  radius,  R25,  at  which  a 
galaxy’s  surface  brightness  falls  to  25  mag  arcsec"^)  of  the  host  galaxies.  Also,  bar 
radii  are  typically  shorter  relative  to  disk  radii  in  late-type  galaxies  than  in  early-type 
galaxies.  In  the  rather  diverse  sample  of  barred  spirals  considered  by  Elmegreen  and 
Elmegreen  (1985),  bars  provide  anywhere  from  -20%  to  -40%  of  the  total  luminosity 
within  the  radius  of  the  bar,  and  from  a few  to  -20%  of  the  total  luminosity  within 
R25.  Moreover,  bars  in  early-type  galaxies  are  generally  stronger,  more  rectangular, 
and  possess  flatter  major  axis  luminosity  profiles  than  bars  in  late-type  galaxies.  The 
latter  tend  to  be  more  elliptical  in  shape  and  possess  light  profiles  that  are  centrally 
peaked  and  fall  off  exponentially  along  the  bar  major  axis. 

Observations  of  the  velocity  fields  of  barred  galaxies  indicate  that  there  exist 
significant  noncircular  streaming  motions  along  the  bar  major  axis  (e.g.  Kormendy 
1983;  Bettoni  et  al.  1988;  Jarvis  et  al.  1988).  This  result  is  consistent  with  the  idea, 
discussed  in  more  detail  below,  that  bar  structure  is  dominated  by  orbits  elongated  along 
the  bar.  Comparison  of  bar  radii  with  the  azimuthally  averaged  rotation  curves  of  their 
host  galaxies  indicates  that  early-type  galaxies  tend  to  have  bars  that  extend  beyond 
the  rising  parts  of  the  rotation  curves,  while  bars  in  late-type  systems  tend  to  end  at  or 
before  the  rotation  curves  peak  or  level  out  (Elmegreen  and  Elmegreen  1985).  Typical 
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values  of  the  central  stellar  velocity  dispersion  are  from  150  to  200  km  s“^,  with  radial 
profiles  of  the  velocity  dispersion  vtu’ying  from  flat  to  sharply  falling  (Sellwood  and 
Wilkinson  1992).  Finally,  estimates  of  the  angular  rotation  rates  of  bars  can  be  made 
by  matching  the  observed  gas  flow  patterns  with  hydrodynamical  models  in  which  the 
bar  tumble  rate  is  varied  (e.g.  Athanassoula  1992). 

Although  much  work  remains  to  be  done  concerning  the  structure  and  evolution 
of  barred  galaxies,  our  understanding  of  the  origin  of  bars  seems  to  be  quite  advanced. 
Some  of  the  earliest  A-body  simulations  of  collisionless,  rotationally  supported  stellar 
disks  dramatically  exemplified  the  bar  instability  (e.g.  Miller  and  Prendergast  1968; 
Hockney  and  Hohl  1969;  Hohl  1971).  Shu  (1970)  and  Kalnajs  (1971)  formulated  the 
problem  of  the  stability  of  axisymmetric  stellar  disks  in  terms  of  normal  modes.  The 
complete  solution,  however,  is  rather  lengthy,  and  the  shapes  and  eigenffequencies  of 
the  unstable  modes  have  been  calculated  only  in  a small  number  of  cases;  nevertheless, 
A-body  simulations  of  a few  of  these  cases  have  yielded  dominant  unstable  modes 
in  close  agreement  with  the  theory  (e.g.  Sellwood  and  Athanassoula  1986).  Toomre 
(1981)  proposed  that  the  bar  instability  is  driven  by  the  positive  feedback  of  swing- 
amplified  waves.  The  idea  here  is  that  any  initial  leading  disturbance  will  propagate  out 
to  corotation,  where  it  will  unwind  and  become  a trailing  disturbance.  As  it  unwinds, 
it  is  swing-amplified  by  the  coordination  of  the  local  epicyclic  frequency  with  the 
rate  of  unwinding  of  the  leading  wave.  The  amplified  trailing  wave  then  propagates 
inward,  and  if  it  is  able  to  reach  the  center,  it  will  be  reflected  as  an  outgoing  leading 
wave,  thereby  starting  the  amplification/reflection  process  over  again.  Wave  action  is 
conserved  at  corotation  by  the  generation  of  a trailing  wave  which  propagates  outward 
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toward  the  outer  Lindblad  resonance  (OLR).  Frequencies  which  cause  the  phases  of  the 
waves  to  match  are  associated  with  standing  waves.  The  interference  of  equal-amplitude 
swing-amplified  leading  and  trailing  standing  waves  leads  to  bar  formation. 

There  are  several  ways  to  control  the  bar  instability.  One,  noted  by  Hohl  (1971), 
is  simply  to  increase  the  velocity  dispersion.  This  tends  not  only  to  stabilize  a 
stellar  disk  against  local  axisymmetric  instabilities  by  raising  the  value  of  Toomre’s 
Q parameter,^  but  also  to  inhibit  the  bar  instability  by  reducing  the  effectiveness  of  the 
swing  amplification  mechanism.  Another  way  is  to  immerse  the  disk  in  a massive  halo. 
This  idea  stems  from  the  work  of  Ostriker  and  Peebles  (1973),  who  found  empirically 
that  their  A-body  models  were  stable  against  bar  formation  if  the  ratio  of  the  rotational 
kinetic  energy  of  the  disk  to  the  total  absolute  potential  energy  was  less  than  a critical 
value  of  0.14  ± 0.02.  A third  way,  suggested  by  Toomre  (1981),  is  to  interrupt  the 
feedback  loop.  This  will  occur  if  the  ingoing  waves  are  prevented  from  reaching  the 
center.  The  existence  of  an  inner  Lindblad  resonance  (ILR),  which  damps  ingoing 
waves  via  wave-particle  interactions,  provides  such  a mechanism.  Toomre  proposed 
that  strong  spiral  perturbations  generated  in  an  unstable  galaxy  shift  material  inward 
and  raise  the  value  of  fl  — f until  an  ILR  appears  and  shuts  off  the  feedback  loop. 

Much  work  has  been  done  concerning  the  structure  and  dynamics  of  bars.  A 

major  portion  of  it  deals  with  the  periodic  orbit  families  that  exist  in  model  bar 

potentials.  Contopoulos  and  Grosbpl  (1989)  give  a comprehensive  review  of  this  topic. 

Sparke  and  Sellwood  (1987)  and  Pfenniger  and  Friedli  (1991)  have  examined  the  orbital 

1.  Toomre  (1964)  found  that  the  condition  for  stability  of  a differentially  rotating  stellar  disk 
to  local  axisymmetric  instabilities  is  (?  = wZqQv  > 1»  where  Or  is  the  local  radial  velocity 
dispersion,  k is  the  local  epicyclic  frequency,  G is  the  Newtonian  gravitation  constant,  and  S 
is  the  local  surface  density. 
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structures  of  numerically  generated  two-dimensional  and  three-dimensional  A^-body  bars, 
respectively.  While  the  work  of  Pfenniger  and  Friedli  presents  some  exciting  new  results 
concerning  the  behavior  of  orbits  and  the  structure  of  bars  perpendicular  to  the  primary 
disk  plane,  we  restrict  our  attention  in  the  present  thesis  to  the  two-dimensional  case. 
The  main  results  of  all  this  work  are  (1)  the  overwhelming  majority  of  particles  that  make 
up  the  bar  are  on  orbits  trapped  around  the  “central”  or,  in  the  notation  of  Contopoulos 
(1975),  Xi  family  (i.e.  the  family  that  reduces  to  circles  in  the  axisymmetric  case)  and 
(2)  bars  end  at,  or  slightly  inside,  corotation.  In  their  study  Sparke  and  Sellwood  found 
a significant  population  of  particles  that  had  sufficient  energy  to  be  able  to  move  freely 
between  the  bar  and  outer  disk  (i.e.  their  Jacobi  constants  exceeded  the  value  at  the 
Lagrange  points  L\  and  L2).  Pfenniger  and  Friedli  also  noted  such  a population  (which 
they  termed  “hot”)  in  their  study. 

The  existence  of  a bar  in  a galaxy  provides  at  least  three  additional  possible 
mechanisms  for  the  generation  of  spiral  structure  (cf.  Elmegreen  and  Elmegreen  1985): 
(1)  a rotating,  growing  bar  may  excite  transient  stellar  spiral  arms  (a  phenomenon 
observed  by  James  and  Sellwood  [1978]  and  Sellwood  [1981]  in  their  numerical 
simulations),  (2)  a rotating,  static  bar  may  excite  spiral  waves  in  the  gas  of  the  outer 
disk  (a  phenomenon  observed,  for  example,  in  the  hydrodynamical  model  of  Sanders 
and  Huntley  [1976]),  and  (3)  a rotating,  static  bar  may  drive  stellar  spirals  in  a strongly 
self-gravitating  disk  (Julian  and  Toomre  1966;  Feldman  and  Lin  1973;  Lin  and  Lau 
1975;  Goldreich  and  Tremaine  1978). 

The  blue  and  near-infrared  surface  photometry  of  fifteen  barred  spiral  galaxies 
provided  by  Elmegreen  and  Elmegreen  (1985)  allows  us  to  begin  to  discriminate 
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between  these  structure-generating  mechanisms.  They  find  that,  in  general,  the  spiral 
arms  of  their  sample  galaxies  contain  star  formation  regions  superposed  on  density 
enhancements  in  the  underlying  stellar  disk.  Also,  they  find  that  the  amplitudes  of 
the  spiral  arms,  as  measured  by  the  near-infrared  brightness  difference  between  the 
arm  and  interarm  regions,  decreases  outward  for  early-type  galaxies  (typically  with 
flat  bars,  as  mentioned  before)  and  increases  outward  for  late-type  galaxies  (which 
generally  have  exponential  bars);  moreover,  they  find  that  galaxies  with  strong  bars 
also  tend  to  have  strong  spirals.  Regarding  spiral  arm  morphology,  they  note  that 
all  of  their  sample  galaxies  that  have  irregular  spiral  arms  also  have  short  bars  of 
exponential  or  indeterminate  type.  On  the  other  hand,  all  galaxies  with  flat  bars  have 
long,  continuous  arms  or  a smooth  ring.  The  Elmegreens  suggest  that  this  apparent 
bimodal  distribution  of  barred  spiral  types  indicates  fundamentally  different  origins  and 
natures.  According  to  this  suggestion,  early-type  barred  spirals  have  bars  that  extend 
close  to  their  corotation  radii  and  continuously  excite  spiral  structure  as  they  slowly 
grow.  Late-types,  on  the  other  hand,  have  rather  weak  bars  that  extend  only  to  the 
ILR  and  excite  correspondingly  weak,  often  multiple,  asymmetric  and  patchy  spiral 
arms.  These  interpretations  have  received  some  recent  support  from  the  numerical 
simulations  of  Combes  and  Elmegreen  (1992). 

Given  the  numerous  possibilities  concerning  the  nature  of  the  spiral  arms  and  the 
positions  of  the  major  resonances  in  barred  spiral  galaxies,  it  is  of  interest  to  ascertain 
whether  or  not  self-consistent  models  of  these  galaxies  can  be  constructed.  Since 
realistic  models  of  galaxies  are,  in  general,  non-integrable,  we  must  resort  to  numerical 
techniques  in  order  to  produce  self-consistent  models.  The  first  effort  of  this  type  was 
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made  by  Martin  Schwarzschild  (1979),  who  succeeded  in  constructing  a self-consistent 
model  of  a nonrotating,  triaxial  stellar  ellipsoid  via  the  method  of  linear  programming. 
His  approach  consisted  basically  of  computing  a library  of  orbits  in  the  potential  derived 
from  an  assumed  density  distribution.  The  time-averaged  density  in  each  cell  of  a grid 
that  spanned  the  volume  of  the  model  was  calculated  for  each  orbit.  The  technique  of 
linear  programming  was  then  employed  to  derive  the  nonnegative  population  weights 
required  in  order  for  the  ensemble  of  orbits  to  give  back  the  imposed  model  density. 
Pfenniger  (1984b)  has  used  a nonnegative  least  squares  technique  to  find  the  population 
weights  for  a self-consistent  model  of  a rapidly  rotating,  two-dimensional  Ferrers  bar. 

The  most  recent  approach  to  the  problem  of  self-consistent  galactic  models  has 
been  developed  by  Contopoulos  and  Grosbpl  (1986,  1988).  Patsis  (1990)  and  Patsis  et 
al.  (1991)  have  successfully  applied  this  method  to  a sample  of  thirteen  unbarred  spiral 
galaxies.  We  extend  that  work  in  this  thesis  by  applying  the  method  of  Contopoulos 
and  Grosbpl  to  a sample  of  three  barred  spirals:  NGC  1398,  NGC  3992,  and  NGC 
1073.  In  Chapter  2 we  describe  in  detail  the  method  of  Contopoulos  and  Grosbpl,  our 
modifications  of  this  method  in  order  to  apply  it  to  the  case  of  barred  spirals,  and  the 
smoothed  particle  hydrodynamics  (SPH)  method  used  to  calculate  the  gas  response.  In 
Chapter  3 we  present  the  relevant  observations  together  with  our  most  successful  models 
of  the  program  galaxies.  In  Chapter  4 we  present  the  results  of  varying  the  parameters 
of  the  most  successful  models.  In  Chapter  5 we  detail  the  role  of  stochastic  orbits  in 
these  models,  and  we  give  the  results  of  the  gas  dynamical  calculations  in  Chapter  6. 
Finally,  we  summarize  the  main  results  in  Chapter  7. 


CHAPTER  2 

MODELING  TECHNIQUES 

The  Method  of  Contopoulos  and  Grosbdl  (CG  Method) 

Before  discussing  the  method  of  Contopoulos  and  Grosbpl  as  applied  to  barred 
spiral  galaxies,  it  is  useful  to  describe  the  technique  as  it  was  originally  applied  to 
unbarred  spiral  galaxies  (Contopoulos  and  Grosbpl  1986,  1988;  Patsis  1990;  Patsis  et 
al.  1991).  The  starting  point  of  any  modeling  attempt  is  the  observational  data.  Here  the 
relevant  data  are  the  axisymmetric  rotation  curve  and  the  surface  photometry,  preferably 
in  the  red  or  near  infrared.  The  former  gives  a rough  estimate  of  the  gravitational  mass 
distribution  in  the  galaxy,  at  least  insofar  as  the  dynamics  of  the  disk  are  concerned. 
The  latter  gives  important  information  not  only  about  the  axisymmetric  distribution  of 
the  stellar  component  of  the  disk,  but  also  about  the  shape  and  strength  of  any  bar 
and/or  spiral  perturbation  which  may  be  present. 

One  kinematic  result  from  observations  of  spiral  galaxies  is  that  the  overwhelming 
majority  of  them  exhibit  rotation  curves  which,  after  steep  initial  rises,  are  flat  or  slowly 
rising  as  far  out  as  they  can  be  determined  (e.g.  Rubin  et  al.  1978;  Bosma  1981).  Any 
reasonable  model  rotation  curve  must  exhibit  this  behavior  as  well.  Secondly,  it  is 
observed  that  some  spiral  galaxies  have  more  prominent  nuclear  bulges  than  others.  In 
fact,  this  feature  has  been  used  as  a classification  criterion  for  these  galaxies  (Sandage 
1961).  Hence  it  is  desirable  to  have  a model  rotation  curve  which  can  represent  the 
bulge  and  disk  contributions  separately.  Lastly,  the  underlying  axisymmetric  potential. 


11 


12 


Vo(r),  derived  from  the  rotation  curve  should  be  expressible  using  standard  analytic 
functions.  The  two  model  rotation  curves  employed  by  Contopoulos  and  Grosbpl  both 
have  these  three  features.  The  one  that  has  been  used  most  extensively  to  date  is 


v(r)  = Umax  exp  (-£fcr)  + [1  - exp(-£rfr)], 


(2-1) 


where  Sb~^  and  ed~^  are  the  scale  lengths  for  the  bulge  and  disk  components,  respec- 
tively. The  quantity  Umax  is  the  asymptotic  value  of  the  circular  rotation  velocity  for 
large  r,  and  the  bulge  fraction  ft  gives  the  importance  of  the  bulge  relative  to  the  disk. 
Figures  2-1  and  2-2  show  the  effects  of  varying  et  and  fb  on  the  shape  of  the  rotation 
curve.  The  axisymmetric  potential  Vo(r),  derived  from  the  relation  is 

given  by 


= -v{^^^{fbexp{-ebr)  - [Inr  + Ei{€dr)]), 


(2-2) 


where  £i(jc)  is  the  first  exponential  integral.  In  this  case  the  angular  velocity  fl(r)  = 
goes  to  infinity  for  r = 0,  thereby  giving  exactly  one  inner  Lindblad  resonance  (ILR).  For 
this  reason  a similar,  but  slightly  modified,  rotation  curve  was  developed  which  allows 
either  zero  or  two  ILRs.  The  circular  rotation  velocity  for  this  model  is  given  by 


u(r)  — Umax 


-Sbv)  -f-  [1  - (1  -I-  £(/r)  exp  (-£(tr)]5 


(2-3) 


while  the  corresponding  axisymmetric  potential  is 


Vo{r)  = -vl^^^{fb{l  + Sbr)exp{-Sbr)  - [Inr  + E\{edr)  + exp(-£rfr)]}.  (2-4) 
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Figure  2-1:  The  effect  of  varying  £b  on  the  shape  of  the  model  rotation  curve  given  by  Eq. 
(2-1):  (1)  £b  = 10e<j;  (2)  £b  = (3)  tb  = '2-£d\  (4)  £b  = £d-  In  all  cases /<,  = 1. 


The  first  step  in  applying  the  CG  method  to  a real  spiral  galaxy,  then,  is  to  fit  one 
of  these  model  rotation  curves  to  the  observed  rotation  curve.  Since  the  observations 
heretofore  have  not  been  detailed  enough  to  allow  one  to  determine  which  model  curve 
is  more  appropriate  for  a given  spiral  galaxy,  the  first  model  has  been  preferentially 
applied  due  to  its  simplicity.  The  fitting  can  be  done  in  any  of  several  ways.  Patsis 
(1990)  has  used  a least  squares  method. 

Once  the  axisymmetric  rotation  curve  v{r)  is  determined  and  the  corresponding 
potential  Vo(r)  derived,  the  next  step  is  to  utilize  the  available  surface  photometry  in 
order  to  find  the  parameters  associated  with  the  spiral  perturbation.  It  is  preferable,  in 
the  present  context,  to  use  red  or  near  infrared  data  because  these  wavelength  ranges 
best  isolate  the  gravitationally  dominant  stellar  components  of  galaxies.  Examples  of 
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the  types  of  parameters  that  can  be  derived  from  the  analysis  of  surface  photometry 
can  be  found  in  Grosbpl  (1985),  Elmegreen  and  Elmegreen  (1985),  Kennicutt  (1981), 
and  Ohta  et  al.  (1990). 


Figure  2-2:  The  effect  of  varying  ft,  on  the  shape  of  the  model  rotation  curve  given  by  Eq. 

(2-1):  (1)  fb  = 0;  (2)fb  = 0.5;  (3)  A = 1;  (4)  fb  = 2.  In  all  cases  eb  = 5sd. 

A cursory  examination  of  the  images  of  spiral  galaxies  indicates  the  general  form 
that  the  spiral  perturbation  must  take.  The  spiral  perturbation  must  increase  from  zero 
outward  from  the  center,  reach  a maximum  somewhere  in  the  disk,  and  then  die  off 
to  zero  again  at  large  distances.  It  must  also  reproduce  the  shape  of  the  arms.  The 
model  perturbation  originally  used  by  Contopoulos  and  Grosbpl  (1986)  to  investigate 
nonlinear  dynamical  effects  at  the  4/1  resonance  has  these  properties.  It  is  a two-armed 
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logarithmic  spiral  perturbation  given  by 


Vi{i\  0)  = /l/'exp  (— £«r)  cos 


(2-5) 


where  i'o  is  the  pitch  angle  of  the  spiral,  is  the  spiral  scale  length,  and  A,  a constant, 
is  a measure  of  the  amplitude.  Also,  the  perturbation  is  assumed  to  be  time-independent 
and  rotating  rigidly  with  a constant  angular  velocity  Qp.  Hence,  implicit  in  this  analysis 
is  the  assumption  of  quasi-stationary  spiral  structure  (QSSS),  a hypothesis  first  put  forth 
explicitly  by  Lin  and  Shu  (1964). 

The  final  form  of  the  spiral  perturbation  that  Contopoulos  and  Grosbpl  (1988)  used, 
and  later  Patsis  (1990)  and  Patsis  et  al.  (1991),  was  altered  by  their  results  of  1986. 
They  found  that  nonlinear  effects  near  the  4/1  resonance  in  galaxies  with  strong  spiral 
perturbations  limit  the  extent  of  any  possible  self-consistent  spiral  structure,  specifically 
to  within  the  4/1  resonance.  Beyond  this  resonance  the  orbits  change  their  orientation 
abruptly  and  the  density  response  maxima  do  not  lie  along  the  imposed  spiral.  This 
necessitated  a change  in  the  form  of  V\{r,0)  to  allow  for  a cutoff  of  arbitrary  abruptness  at 
a specified  radius.  Also,  the  formula  was  generalized  to  allow  for  40  and  60  components 
in  addition  to  the  20  component.  The  final  form  of  the  imposed  spiral  perturbation  can 
be  written  most  generally  as 


Vi{r,0)=  ^ A„,(r)r  exp  (-£3, „r)  cos 


m In  r 


m=2,4,6 


tan  ioi 


— mO  -f  mOfn  I , (2-6) 


m 


where 


1 


A/7i(?')  — s (A,n  A,-,7))_^(1  -t*  tanh  [kojk  (* '2  m ^)])  “b  -^Tm 


(2-7) 


and  Am,  Arm,  K2m  and  r2m  are  all  constants.  Am  measures  the  maximum  amplitude 
achieved  by  the  mO  component.  Arm  gives  any  residual  amplitude  of  the  mO  component 


16 


that  extends  beyond  its  cutoff  radius  r2„„  and  K2m  defines  the  steepness  of  the  cutoff. 
By  analogy  with  equation  (5),  esm~^  is  the  scale  length  and  /omthe  pitch  angle  of  the 
mO  component.  Also,  differences  in  the  phases  of  the  various  components  are  allowed 
with  the  corresponding  values  of  Om- 

Analysis  of  the  surface  photometry  can  yield  (at  least  rough)  estimates  of  most 
or  all  of  the  parameters  needed  to  specify  the  spiral  perturbation.  It  is  not  necessary 
at  this  stage  to  find  the  specific  values  of  the  parameters  which  yield  the  most  self- 
consistent  model.  Variations  of  and  corrections  to  the  parameters  are  made  after  the 
results  of  subsequent  modeling  attempts  are  known.  In  fact,  this  is  the  essence  of  the 
CG  method.  The  observations  do,  however,  place  limits  on  the  ranges  of  acceptable 
parameter  values,  as  well  as  provide  good  initial  guesses. 

After  the  parameters  of  the  spiral  perturbation  are  determined,  at  least  roughly, 
from  the  surface  photometry,  the  next  step  is  to  calculate  the  periodic  orbits.  For  this 
the  equations  of  motion  are  needed.  Since  these  models  are  restricted  to  two  dimensions 
(i.e.  the  disk  plane),  the  equations  of  motion  can  be  derived  from  the  two-dimensional 
Hamiltonian 

//  = i - n„Jo  + v^r)  + Viir,  0)  = h,  (2-8) 

where  r is  the  radial  component  of  the  velocity,  Jq  the  angular  momentum  in  the  inertial 
frame  of  reference,  Q,p  the  angular  velocity  of  the  reference  coordinate  system  (i.e.  the 
frame  corotating  with  the  spiral  pattern)  with  respect  to  the  inertial  frame,  and  h the 
numerical  value  of  H.  Also,  Vo(r)  and  V\{r,0)  are  the  axisymmetric  and  perturbation 
potentials,  respectively.  From  Eq.  (2-8)  the  necessary  equations  of  motion  can  be 
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derived.  They  are,  in  polar  coordinates  in  the  rotating  frame. 


M 

dt 


j],), 


dr 

dt 


(2-9) 


These  equations  can  be  solved  numerically  by  any  of  the  standard  integration  techniques. 
The  one  that  has  been  used  by  Contopoulos  and  Grosbpl  (1986,  1988)  is  a fourth-order 


Runge-Kutta  method  with  a variable  stepsize,  which  allows  the  relative  error  in  each 
variable  to  be  maintained  below  a preset  level  (typically  one  part  in  10^).  The  error 
controls  allow  for  energy  conservation  along  the  integrated  orbit  to  the  order  of  one 
part  in  10^  for  a typical  orbit  integration. 

The  periodic  orbits  of  particular  interest  here  are  the  ones  comprising  the  “central,” 
or  in  the  notation  of  Contopoulos  (1975),  Xi  family,  which  reduce  to  circular  orbits  in 
the  purely  axisymmetric  case.  The  introduction  of  a bar  or  spiral  perturbation  breaks 
this  central  family  into  an  infinity  of  families  by  gaps  at  all  even  resonances  between  the 
epicyclic  frequency  «:(r),  defined  by  k“{v)  = - - -I- 1 = 4fl'^(r)  + 2rfl(r)^^^, 

and  the  orbital  angular  frequency  measured  in  the  frame  corotating  with  the  perturbation 
Q(r)  - flp  (Contopoulos  1983).  Here  H(r)  = That  is,  where 


fl(r)  - Qj, 


(2-10) 


Also,  regions  of  instability  are  produced  in  the  central  family  near  the  odd  resonances, 
where 


/v(r) 

H(?')  — flp 


(2-11) 


From  these  unstable  regions  bifurcate  the  odd  resonant  families. 
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Since  the  available  dynamical  evidence  (Contopoulos  and  Grosb0l  1986,  1988; 
Patsis  1990;  Patsis  et  al.  1991)  indicates  that  the  spiral  structure  of  normal  spiral  galaxies 
extends  from  the  ILR  (2/1  resonance)  to  either  the  4/1  resonance  (for  strong  spirals)  or 
the  corotation  resonance  (for  weak  spirals),  one  can  get  a reasonable  estimate  of  the 
pattern  speed  Qp.  This  parameter  is  notoriously  difficult  to  determine  observationally, 
and  therefore  must  be  left  adjustable  in  the  fitting  procedure.  Nevertheless,  observations 
of  the  extent  of  the  spiral  structure,  together  with  the  axisymmetric  rotation  curve,  can 
produce  reasonable  initial  values. 

The  significant  branches  of  the  central  family,  then,  that  exist  in  the  fully  perturbed 
model  between  the  inner  Lindblad  resonance  and  corotation  are  found.  In  this  region 
there  is  a one-to-one  correspondence  between  the  radius  rc  of  the  circular  periodic  orbit 
of  the  axisymmetric  case  and  its  value  of  the  Hamiltonian  h.  This  allows  for  a rather 
simple  algorithm  to  find  the  periodic  families  of  the  fully  perturbed  case.  For  example, 
consider  the  branch  of  the  central  family  which  exists  between  the  ILR  and  the  4/1 
resonance.  The  radius  of  the  periodic  orbit  halfway  between  these  resonances  in  the 
axisymmetric  case  will  have  the  value  ro  = rc  = 7(riLR  + ^4/1)5  where  riLR  and  r4/i 
are  the  positions  of  the  inner  Lindblad  and  4/1  resonances,  respectively.  The  orbit  will 
cross  the  0 = 0 axis  with  a radial  velocity  kq  = 0.  Now  let  the  “amplitude”  of  the 
perturbation  potential  be  increased  from  zero  by  some  small  fraction  of  its  maximum 
value.  Also,  let  the  value  of  the  Hamiltonian  for  the  new  orbit  remain  h = h{rc).  If 
this  new  orbit  is  started  at  radius  ro  with  zero  radial  velocity,  however,  it  will  not,  after 
being  integrated  for  one  revolution,  cross  the  = 0 axis  at  radius  ro  with  zero  radial 
velocity.  That  is,  it  will  not  be  periodic  due  to  the  imposed  perturbation.  But  since  the 
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applied  perturbation  is  small,  the  orbit  will  be  close  to  a periodic  orbit,  thus  the  actual 
periodic  orbit  can  be  readily  determined  using  a two-dimensional  Newton  method. 

Consider  an  arbitrary  orbit  of  fixed  “energy”  h,  close  to  a periodic  orbit,  starting 
with  initial  conditions  (ro,  ro)  along  the  0 = 0 axis.  After  one  revolution  it  will  cross 
the  9 = 0 axis  again  but  with  different  radius  and  radial  velocity  (ri,  ri).  Now  vary 
separately  ro  and  ro  by  Ar  and  Ar,  respectively,  and  integrate  these  new  orbits  for  one 
revolution  each,  producing  the  following  mappings  of  initial  to  final  conditions:  (ro  -i- 
Ar,  ro)  — ^ (r2,  K2)  and  (ro,  ro  + Ar)  (r^,  rs).  These  quantities  are  now  used  to 
compute  the  partial  derivatives  of  the  final  radius  and  radial  velocity,  Vf  and  r/,  with 
respect  to  the  initial  radius  and  radial  velocity,  r/  and  r,-.  Specifically, 


dvf  V2  — ri 

drf  T3  — r\ 

dri  Ar  ’ 

dri  Af  ’ 

drf  r2  — ri 

drf  rs  — r\ 

(2-12) 


dvi  Ar  ’ dri  Ar 

These  numerically  determined  partial  derivatives  are  used  to  find  the  first-order  correc- 
tions to  ro  and  ro  necessary  to  adjust  ro  and  ro  more  closely  to  the  initial  conditions  of 
the  periodic  orbit.  The  corrections,  8r  and  8r,  are  determined  by  solving  the  following 


equations: 


dvf  drf 

n + -::;—8r  -f  -Trr-or  = ?’o  + 8r 
on  on 


(2-13) 


drf  dvf 

fi  -f  -7^8r  -f  ~^8r  - ?'o  -f  8r, 

UT  i UT I 


These  corrections  are  added  to  ro  and  ro. 


ro  — > ro  -I-  8r 
ro  — ^ ro  + 8r, 


(2-14) 
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thereby  improving  the  estimates  of  the  initial  conditions  of  the  periodic  orbit.  This 
procedure  is  then  repeated  until  |^r|  and  l^fj  fall  below  some  predefined  tolerance 
level.  There  is  usually  no  problem  with  convergence  of  this  method  if  the  starting 
initial  conditions  are  not  too  “far”  from  the  periodic  orbit. 

In  this  manner  it  is  possible  to  “step  up”  the  perturbation  amplitude  from  zero  to 
its  full  value  in  small  increments,  finding  accurately  (typically  to  one  part  in  10^®)  the 
periodic  orbit  at  each  step.  After  the  first  step,  better  starting  guesses  at  the  initial 
conditions  of  the  periodic  orbit  of  the  next  step  are  provided  by  an  extrapolation 
procedure.  This  whole  process  is  done  for  only  one  orbit  between  any  given  pair 
of  even  resonances  inside  corotation.  The  orbit  “energy”  h,  parameterized  by  rc,  is 
chosen  to  be  that  of  the  circular  orbit  (in  the  axisymmetric  case)  halfway  between  the 
positions  of  the  resonances  in  question,  as  in  the  example  of  the  ILR  and  4/1  resonance 
above.  Once  this  periodic  orbit  is  found  in  the  fully  perturbed  case,  the  whole  family 
of  orbits  between  the  resonances  is  calculated  by  varying  the  energy  rc  both  positively 
and  negatively.  Again,  after  the  first  step  each  way,  an  extrapolation  procedure  is  used 
to  provide  better  first  guesses  at  the  initial  conditions  of  the  periodic  orbit. 

In  models  of  unbarred  spiral  galaxies,  the  resonance  ranges  probed  for  branches  of 
the  x\  family  are  usually  limited  to:  2/1  (ILR)  to  4/1,  4/1  to  6/1,  and  6/1  to  8/1.  This  is 
because  in  a realistic  model  of  a spiral  galaxy,  the  positions  of  the  resonances  bunch  up 
near  corotation  and  most  of  the  galaxy  inside  corotation  is  also  inside  the  8/1  resonance. 
For  example,  in  the  standard  model  for  NGC  5247  considered  by  Contopoulos  and 
Grosbpl  (1986),  corotation  is  roughly  at  23  kpc.  In  comparison,  the  ILR  is  around  1.5 
kpc,  the  4/1  resonance  around  12  kpc,  the  6/1  resonance  around  16  kpc,  and  the  8/1 
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resonance  around  17.5  kpc.  And,  according  to  Contopoulos  and  Grosbpl  (1986,  p.l55), 
“the  congestion  of  infinite  resonances  n/1  near  corotation  produces  a large  degree  of 
stochasticity  there.”  Still,  in  cases  where  significant  stable  periodic  families  exist  beyond 
the  8/1  and/or  10/1  resonances,  orbits  in  these  regions  are  generally  included. 

Now  that  the  families  of  periodic  orbits  to  be  included  have  been  determined, 
the  next  step,  undoubtedly  the  most  important  step  in  the  method  of  Contopoulos  and 
Grosbpl,  is  to  generate  the  density  response  map.  For  this  a two-dimensional  polar 
grid  is  employed;  moreover,  since  the  entire  problem  is  symmetric,  only  half  of  the 
grid  is  actually  used  in  order  to  save  computational  time  and  storage  space.  Use  of  a 
polar  grid  permits  both  a simple  way  of  analyzing  the  self-consistency  of  the  density 
response,  and  higher  resolution  near  the  center  where  the  density  and  density  gradient 
become  large.  The  radial  range  of  the  grid,  as  well  as  the  radial  and  azimuthal  widths 
of  the  grid  cells,  can  be  varied  in  the  analysis.  The  only  restriction  is  that  the  number 
of  azimuthal  cells  should  be  a power  of  2.  The  reason  for  this  will  be  outlined  later  in 
the  discussion  of  the  analysis  of  the  density  response. 

The  density  response  map,  then,  is  constructed  by  computing  a representative 
sample  of  the  periodic  orbits,  and  a set  of  orbits  dispersed  in  radial  velocity  around 
them,  and  incrementing  the  value  of  the  density  of  a given  cell  by  the  product  of  the 
weight  of  the  orbit,  to  be  discussed  in  detail  shortly,  and  the  amount  of  time  spent  by 
the  orbit  in  the  cell.  Each  periodic  orbit  family  included  in  the  analysis  is  sampled 
at  fixed  intervals  of  “energy,”  as  parameterized  by  the  radius  of  the  circular  orbit 
of  equivalent  energy  in  the  axisymmetric  case.  The  sampling  step  in  the  “energy”  so 
parameterized  is  taken  to  be  times  the  radial  width  of  a grid  cell,  where  NOR  = 
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1,  2,  3,  ...  is  an  input  parameter.  In  all  of  the  models  considered  in  this  thesis  we  have 
taken  the  value  of  NOR  to  be  2. 

Every  periodic  orbit  chosen  by  the  sampling  process  to  be  included  in  the  density 
map  is  then  integrated  for  one  half  period  (since  only  half  of  the  polar  grid  is  used)  in 
order  to  have  its  contribution  to  the  density  map  added  and  to  generate  initial  conditions 
for  the  dispersed  orbits.  The  periodic  orbits,  as  originally  determined,  all  have  their 
initial  conditions  (r,  r)  specified  along  the  0 = 0 axis.  It  is  desirable,  however,  in  order 
to  estimate  more  accurately,  and  with  reduced  integration  times  per  orbit,  the  actual 
density  response  of  orbits  dispersed  around  the  periodic  orbit,  to  start  orbits  for  the 
density  response  calculation  away  from  this  axis  but  along  the  trajectory  of  the  periodic 
orbit.  Therefore,  during  the  integration  of  the  periodic  orbit,  its  radius  r,  radial  velocity 
r,  azimuth  6,  and  angular  momentum  Jq  are  stored  every  fraction  of  its  half  period, 
where  NOA  = 1,  2,  3,  ...  is  also  an  input  parameter. 

After  the  density  contribution  of  each  periodic  orbit  is  added  to  the  density  map, 
the  orbits  dispersed  about  the  periodic  orbit  are  computed.  Whereas  the  periodic  orbit 
is  calculated  for  its  full  half  period,  it  is  not  really  necessary  to  calculate  the  dispersed 
orbits  for  this  length  of  time  due  to  the  increased  numbers  of  dispersed  orbits  calculated 
per  periodic  orbit.  The  actual  time  of  integration  for  the  dispersed  orbits  is  controlled 
by  a parameter  PN  such  that 


integration  time  = to  = PN  x 


1 PZ 
NOA  ^ 


(2-15) 


where  PZ  is  the  period  of  the  base  periodic  orbit.  The  value  of  PN  is  typically  chosen 
to  be  of  order  NOA,  empirically  determined  to  be  a good  compromise  between  the  need 
to  generate  an  adequate  representative  density  response  for  the  orbit  and  the  need  to 
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minimize  computational  time  per  orbit.  The  dispersed  orbits  are  started  at  the  NOA 
stored  points  along  the  periodic  orbit.  The  distribution  of  dispersed  orbits  is  assumed  to 
be  Gaussian  with  velocity  dispersion  a,  thus  the  orbits  must  be  dispersed  and  weighted 
appropriately  relative  to  the  central  periodic  orbit  in  order  to  estimate  best,  using  a finite 
number  of  orbits,  a continuous  Gaussian  distribution.  A quadrature  scheme  is  indicated. 

Consider  a periodic  orbit,  characterized  by  (r,  9,  r,  Jq),  a point  along  the  trajectory 
of  the  orbit.  The  problem  is  to  compute  the  surface  response  density  of  a continuum  of 
orbits  all  starting  at  this  point,  but  distributed  in  radial  velocity  about  r in  a Gaussian 
fashion  with  dispersion  a.  Let  the  surface  response  density  of  the  orbit  started  at  (r,  6, 
r,  Jo)  and  calculated  for  time  ro  be  given  by  s{r,6,r,  Jo^to).  The  problem,  then,  can 
be  stated  mathematically  as 


OO 

S{r,6,Jo,to)—  J exp 

— OO 


Jo,  to) dr  , 


(2-16) 


where  S{r,  6,  Jo,  to)  represents  the  total  surface  response  density  of  all  orbits  distributed 
about  the  period  orbit  and  calculated  for  time  to-  In  order  to  estimate  the  above  integral, 
a quadrature  scheme  employing  the  Hermite  polynomial  Hn(x)  is  used  (in  all  cases 
considered  in  this  thesis  we  have  taken  n = 9).  Specifically  (Abramowitz  and  Stegun 


1965,  p.  890), 

J n 

/ exp  (-x^)f{x)dx  = ^ wvif(xi)  + Rn, 

-OO 

where  xi  is  the  zero  of  Hnix),  the  weights  wvi  are  given  by 

'n?[Hn-i{xi)f‘'‘ 


(2-17) 


(2-18) 


and  Rn  is  a small  remainder. 
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Since  we  use  a nine  point  (n  = 9)  Hermite  integration  to  calculate  the  effect  of 
velocity  dispersion  on  the  density  response,  a total  of  nine  orbits,  the  one  periodic  orbit 
and  eight  dispersed  about  it,  four  on  either  side  in  r space,  are  integrated,  with  the 

OO 

initial  conditions  given  in  Table  2—1.  Since  f exp  (^—x^)dx  = y/w,  the  weights  wvi 

— OO 

must  be  normalized  by  dividing  by  a/x  before  the  integrations  are  done. 


Table  2-1:  Initial  conditions  of  the  periodic  orbit  and  the  eight 
nonperiodic  orbits  used  to  calculate  the  effect  of  the  velocity  dispersion. 


Orbit  Number 

Initial  Conditions 

Orbit  Weight 

1 (periodic  orbit) 

r,  6>,r  + (jxi,  Jo 

WVi 

2,3 

r,  0,  r ± ax2,  Jo 

WV2 

4,5 

r,0,r±  ax3,Jo 

wvs 

6,7 

r,  0,  r ± (7X4,  Jo 

WV4 

8,9 

r,0,f±  (7X5,  Jo 

WVs 

Table  2-2:  The  values  of  the  xi  and  wv{ 

of  Table  2-1. 

• 

1 

Xi 

WVi 

1 

0.00000 

7.20235... 

X 

10-1 

2 

0.72355... 

4.32651... 

X 

10-1 

3 

1.46855... 

8.84745... 

X 

10-2 

4 

2.26658... 

4.94362... 

X 

10-2 

5 

3.19099... 

3.96069... 

X 

10-2 

In  addition  to  the  Hermitian  weight  mentioned  above,  we  must  also  include  some 
further  weighting  factors  to  measure  properly  the  contribution  of  a particular  orbit  to 
the  response  density  map.  These  weights  are  introduced  to  account  for  (1)  the  relative 
abilities  of  different  periodic  families  coexisting  at  “energy”  rc  to  trap  matter  around 
them,  (2)  the  population  of  a particular  orbit  of  “energy”  Vc  according  to  the  imposed 
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axisymmetric  surface  density  at  radius  Vc,  (3)  the  “grid”  effects  of  the  parameters  NOA 
and  NOR,  and  (4)  the  time  of  orbit  integration. 

Sometimes,  in  the  course  of  calculating  the  periodic  orbit  families  to  include  in 
the  model,  we  find  more  than  one  stable  family  at  the  same  “energy”  Therefore  the 
question  arises  as  to  the  relative  importance  of  each  orbit  family  in  the  energy  intervals 
where  there  is  overlap.  Contopoulos  (1979)  has  given  an  estimate  of  the  relative  amount 
of  trapping  done  by  a periodic  orbit  of  a given  energy  r^.  This  estimate  is  measured, 
for  the  m*  orbit  family,  by  the  weight 


K“{rc)x 


2 < f2  > 


(2-19) 


where  n{rc)  is  the  epicyclic  frequency  at  radius  Vc,  Xm  is  the  maximum  deviation  of  the 
periodic  orbit  from  r<;,  and  <p->  is  the  square  of  the  RMS  radial  velocity  dispersion. 
Since  the  velocity  dispersion  a mentioned  earlier  is  assumed  to  be  isotropic,  <J^>  = 
Therefore  the  relative  weight  Wm  of  the  m*  family  at  “energy”  rc  is  given  by 


(2-20) 


where  n is  the  total  number  of  stable  orbit  families  coexisting  at  energy  rc.  In  energy 
intervals  where  only  a single  orbit  family  exists,  this  weighting  factor  is  obviously 
equal  to  one. 

The  next  weighting  factor  to  be  considered  is  the  most  obvious.  Orbits  of  lower 
energy  rc,  which  exist  nearer  to  the  center  of  the  galaxy,  should  be  weighted  more  than 
orbits  of  higher  energy,  which  exist  in  the  outer  regions  of  the  galaxy.  This  simply 
reflects  the  fact  that  the  matter  density  in  the  disk  decreases  outward.  Therefore  a 
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weight  wv,  of  the  form 


• wi:  = S(rc),  (2-21) 

is  used.  Here  S(rc)  is  the  surface  density  of  the  unperturbed  disk  at  radius  rc.  This 
quantity  cannot  simply  be  derived  from  the  axisymmetric  rotation  curve,  because  it  is 
not  clear  that  the  majority  of  the  gravitational  mass  in  spiral  galaxies  is  concentrated  in 
or  near  the  observed  disk.  In  fact,  observed  disk  luminosity  profiles  typically  decrease 
exponentially  outward  (Freeman  1970).  Unless  the  actual  mass-to-light  ratio  varies 
wildly  across  the  disk,  it  is  plausible  to  infer  that  the  matter  density  in  the  disk 
also  exhibits  an  exponential  decrease  with  radius,  and  that  the  excess  gravitational 
mass  responsible  for  the  observed  rotation  curve  is  distributed  in  a spherical  or  oblate 
spheriodal  distribution,  but  away  from  the  disk  in  any  case.  Given  these  facts,  the 
surface  density  of  the  unperturbed  disk  is  taken  to  be  exponential,  given  by 


S(r)  = coexp(-£or),  (2-22) 

where  cq  is  the  central  surface  density  and  the  scale  length  of  the  disk.  Thus  the 
weight  factor  applied  for  an  orbit  of  energy  is 

wj:  = Co  exp  (-£o^c)-  (2-23) 


The  final  weighting  factors  are  introduced  in  order  to  eliminate  the  effects  of  the 
particular  choices  of  NOR,  the  number  of  sampling  steps  in  “energy”  Kc  per  grid  cell, 
NOA,  the  number  of  different  starting  positions  for  dispersed  orbits  along  the  periodic 
orbit,  and  to,  the  time  of  integration.  The  grid  effects  are  taken  out  by  a weight  of 
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the  form 


1 


1 


tOgj'id 


X 


NOR  NOA 


(2-24) 


Also,  since  the  contribution  of  a particular  orbit  to  any  one  grid  cell  is  the  previous 
weight  of  the  orbit  multiplied  by  the  amount  of  time  spent  by  the  orbit  in  the  cell,  the 
orbit  weight  must  be  divided  by  the  total  time  of  orbit  integration.  As  noted  before,  the 
integration  time  is  ^PZ  for  periodic  orbits  and  ^PZ  x PN  x j;;j^  for  orbits  dispersed 
about  the  periodic  orbit.  Again,  PZ  is  the  period  of  the  base  periodic  orbit  and  PN 
the  parameter  which  controls  the  length  of  integration  for  nonperiodic  orbits  through 
Eq.  (2-15).  Therefore  the  corresponding  weights  applied  are  simply  the  inverses  of 


these  times. 


^time 


2 

PZ’ 
2xNOA 


periodic  orbits 


(2-25) 


M X P Z ’ nonperiodic  orbits 


The  total  weight  applied  to  a given  orbit,  then,  can  be  written  as 


W^tot  = WV  X Wm  X Wj:  X tCgrid  X tCtime 


(2-26) 


The  complete  density  response  map,  then,  is  produced  by  (1)  specifying  the 
minimum  (RMIN)  and  maximum  (RMAX)  radial  extent  of  the  grid,  as  well  as  the 
radial  (DRS)  and  azimuthal  cell  widths,  (2)  specifying  the  parameters  NOR,  NOA,  and 
PN  described  above,  (3)  specifying  the  central  value,  ao,  and  slope,  ar,  of  the  velocity 
dispersion  (the  form  of  the  velocity  dispersion  is  assumed  to  be  a(r)  = ao  + rur,  a 
linear  relation  being  the  simplest  functional  form  possible  which  allows  a{r)  to  vary 
across  the  disk),  and  (4)  stepping  in  radius  Vc  from  RMIN  to  RMAX  in  increments  of 
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and  integrating  all  orbits  existing  at  each  value  of  rc,  weighted  as  described  above. 
The  results  from  this  raw  map  are  then  preprocessed  slightly  before  being  stored  for 
later  analysis.  Each  semiannulus  of  the  grid  is  processed  in  turn.  Thus,  for  a fixed 
semiannulus,  the  density  value  of  a each  azimuthal  cell  is  “smoothed”  with  the  values 
of  the  two  adjacent  cells  according  to  the  smoothing  formula 

DEN(i)  = iDEN(i)  + i(DEN(i  - 1)  + DEN(i+l)),  (2-27) 

2 d 


where  i is  the  azimuthal  cell  index.  At  the  same  time,  the  total  “density,”  RHO,  in 
the  semiannulus  is  computed. 


n 


RHO  = J]DEN(i), 


(2-28) 


where  n is  the  number  of  azimuthal  cells.  The  value  of  each  cell  is  then  normalized  by 
the  “average”  cell  value  of  the  semiannulus  and  finally  decreased  by  one,  such  that 


nxDEN(0  , 


(2-29) 


represents  directly  the  fractional  deviation  of  the  azimuthal  cell  of  the  semiannulus 
from  the  average. 

The  final  step  of  the  method  of  Contopoulos  and  Grosbpl  is  to  analyze  the 
density  response  for  consistency  with  the  imposed  density.  In  other  words,  the  self- 
consistency  of  the  model  is  measured.  Here  again,  the  analysis  is  done  on  the  semiannuli 
separately.  Two  basic  quantities  are  used  to  estimate  the  self-consistency  of  the  model. 
The  first  compares  the  amplitude  of  the  29  component  of  the  response  density  to  the 
corresponding  amplitude  of  the  imposed  density.  The  second  measures  the  difference 
between  the  phases  of  the  response  and  imposed  density  maxima  in  the  semiannulus. 
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First,  a Fast  Fourier  Transform  (FFT)  with  respect  to  azimuth  is  performed  on  the 
equally  spaced  data  in  array  DEN.  The  reason  for  requiring  the  number  of  azimuthal 
cells  to  be  an  integer  power  of  2 now  becomes  clear.  While  the  FFT  can  be  made  to 
accomodate  data  sets  of  other  lengths,  the  algorithm  is  applied  most  easily  to  a set  of 
length  n = 2*,  yt  a nonnegative  integer.  The  power  and  phase  of  each  positive  frequency 
component  fi  = ^,  / = 0,  1,  2, ...,  nil,  are  computed.  For  the  i = 1 frequency  component 
(i.e.  the  26  component),  the  difference  A0  between  the  response  density  maximum  and 
imposed  potential  minimum  is  then  calculated.  This  angle  difference  is  measured  with 
respect  to  the  imposed  potential  minimum  rather  than  the  imposed  density  maximum 
because  these  two  coincide  almost  exactly  and  because  the  angle  of  the  the  imposed 
potential  minimum  is  more  readily  ascertained.  In  a fully  self-consistent  model, 
should  be  zero  for  all  values  of  the  radius.  Therefore  A^  constitutes  one  measure  of 
the  self-consistency  of  a given  model. 

The  amplitude  of  the  / = 1 frequency  component  corresponds  to  the  quantity 


^2,resp 

^0,resp 


(2-30) 


where  a2,resp  is  the  actual  amplitude  of  the  26  component  of  the  response  density  in 
the  semiannulus  under  consideration  and  ao^resp  is  the  mean  surface  response  density 
in  the  semiannulus,  or  RHO/n.  It  is  of  interest  to  compare  this  with  the  corresponding 
quantity  of  the  imposed  density,  or 


^2, imp 
^0,imp 


(2-31) 


The  mean  imposed  density  at  the  radius  of  the  semiannulus  is  given  by  the  formula  for 
the  exponential  disk  [Eq.  (2-23)].  The  amplitude  of  the  26  component  of  the  imposed 
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density  is  derived  from  the  imposed  potential  via  Poisson’s  equation.  To  this  end 

Vandervoort  (1970;  see  also  Contopoulos  and  Grosbpl  1988)  has  given  the  following 

formulae,  valid  in  the  case  of  a tightly  wrapped  spiral  density  wave,  relating  the  full 

three-dimensional  potential  Vi  and  density  p\  of  the  spiral: 

= yiosechl^'^l(^//A), 


(2-32) 


+ 0(i), 

where  Vio  is  the  spiral  potential  formula  in  the  disk  plane  (z  = 0),  A is  the  z thickness  of 
the  spiral  and  k is  the  wave  number.  The  parameter  A,  then,  allows  the  vertical  thickness 
of  the  spiral  to  be  controlled  and  must  be  included  in  the  adjustment  of  parameters.  For 
the  form  of  spiral  potential  used  by  Contopoulos  and  Grosbpl, 

k = — (2-33) 

r tan  ^o 

where  m is  the  number  of  spiral  arms.  From  the  above  relations  (J2,imp  can  be  computed. 
Specifically, 

OO  OO 

<72, imp  ~ J Pi{z)dz  = 2 J pi{z)dz.  (2-34) 


— OO  0 

Here  A is  taken  to  be  constant  over  the  entire  disk  of  the  galaxy, 
quantity 


<T2 


res] 


(^0,resp 


CT2. 


tmi 


*7’0,t  mp 


At  this  point  the 


(2-35) 


can  be  formed.  In  a fully  self-consistent  model  R*  should  be,  in  principle,  equal  to 
unity  for  all  values  of  the  radius.  In  practice,  however,  it  will  not  be  so  because  ao,resp 
and  (TQ^imp  are  calculated  in  the  different  ways  described  above.  The  actual  value  of  R 
can  be  scaled  over  a rather  wide  range  by  varying  cq,  the  central  surface  density  of  the 
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“imposed”  disk.  The  maximum  value  of  cq  is  limited  only  by  the  total  axisymmetric 
potential.  Therefore,  practical  self-consistency  requires  that  R simply  be  constant  with 
radius.  This  is  the  second  measure  of  the  self-consistency  of  a given  model. 

The  whole  procedure  outlined  above  is  then  repeated  for  different  sets  of  model 
parameters,  and  a “best”  model  achieved  by  a sort  of  manual  iteration.  It  would  be 
desireable  to  automate  the  whole  procedure,  but  at  this  point  too  much  human  decision- 
making is  required  at  intermediate  steps  for  this  to  be  feasible. 

The  CG  Method  Modified  for  the  Case  of  Barred  Spiral  Galaxies 


In  this  section  are  described  the  modifications  made  to  the  CG  method  in  order 
to  apply  it  to  the  case  of  barred  spiral  galaxies.  Perhaps  the  most  significant  change 
is  the  nature  of  the  assumed  model.  In  contrast  to  the  original  CG  method  where  a 
simple  function  is  fit  to  the  observed  axisymmetric  rotation  curve,  thereby  allowing 
the  total  axisymmetric  potential  to  be  derived  in  a correspondingly  simple  form,  model 
components  are  used.  This  is  necessitated  by  the  fact  that  the  full  complement  of 
nonaxisymmetric  structure  cannot,  in  this  case,  be  adequately  modeled  as  a simple 
perturbation  of  the  axisymmetric  disk.  This  approach  was  tried  initially,  by  replacing  the 
quantity  -,-■-■1^-- — m0-\-m9m  with  md  in  the  formula  for  V\{r,9)  for  values  of  radius  less 
than  the  length  of  the  bar,  with  somewhat  limited  success.  It  was  found  that  even  with 
the  inclusion  of  significant  higher  harmonics  the  imposed  bar  density  distribution  could 
not  adequately  simulate  realistic  bars,  even  qualitatively.  Therefore  a three  component 
model  is  employed,  with  components  introduced  to  represent  separately  a bar,  a disk. 


and  a halo. 
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The  axisymmetric  component  is  modeled  as  the  superposition  of  disk  and  halo 
components.  In  view  of  the  fact,  mentioned  above,  that  the  observed  brightness 
profiles  of  disk  galaxies  are  typically  exponential,  an  exponential  disk  is  used.  The 
corresponding  potential  in  the  disk  plane  is  readily  derived  (cf.  Binney  and  Tremaine 
1987,  p.77)  and  can  be  written 

Voir)  = -7cGcQr\lQ{y)Ki{y)  - Ii{y)Ko{y)],  (2-36) 


where  G is  the  Newtonian  gravitation  constant,  Co  is  the  central  surface  density  of 
the  disk,  y = \ear,  and  the  /„  and  Kn  are  modified  Bessel  functions  of  the  first  and 
second  kind,  respectively.  Again,  eo~^  is  the  scale  length  of  the  exponential  disk.  The 
circular  rotation  velocity  V£>{r)  corresponding  to  this  potential  is  readily  derived  from 
the  relation 


dVp{r) 

dr 


(2-37) 


using  the  fact  (Abramowitz  and  Stegun  1965,  p.376)  that 


dlp{x) 

dx 


= h(x) 


dKo(x) 

dx 


= -Ki{x), 


dh{x)  dKi 

= Io[x)  , and 


dx 


dx 


1 

—Ko{x) Ki{x). 

X 


(2-38) 


The  result  is 


Vpir)  = yj i:GcQeQr‘^[h{y)KQ{y)  - h{y)K\{y)].  (2-39) 

The  general  form  of  this  rotation  curve,  along  with  the  Keplerian  curve  for  a point  mass 
equal  to  the  mass  Md  = of  the  exponential  disk  is  shown  in  Figure  2-3  (essentially 
Figure  2-17  of  Binney  and  Tremaine  1987,  p.  78). 
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Figure  2-3:  Rotation  curves  of  the  exponential  disk  (solid 
curve)  and  a point  with  the  same  total  mass  (dotted  curve). 


As  can  be  readily  seen  from  Figure  2-3,  the  rotation  curve  of  an  exponential  disk 
declines  after  approximately  two  disk  scale  lengths  and  rather  quickly  approaches  the 
Keplerian  value.  This,  as  described  above,  is  in  contradiction  to  the  observed  rotation 
curves  of  most  disk  galaxies.  Hence  a separate  halo  component  is  included  in  order 
to  “hold  up”  the  rotation  curve  at  large  radii.  The  density  distribution  chosen  for  this 
component  is  that  of  a Plummer  sphere: 


pjj{r)  = 


(2-40) 


Here  M//  is  the  total  halo  mass  and  b is  a parameter  which  controls  the  central 
condensation.  The  potential  corresponding  to  this  density  distribution  is  particularly 
simple: 

GMh 

\/r^  + 6“ 


(2-41) 
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From  the  potential  the  circular  rotation  velocity  is  derived: 


(2-42) 


Figure  2-4  shows  the  rotation  curve  of  the  Plummer  sphere  and  the  corresponding 
Keplerian  rotation  curve  for  a point  mass  equal  to  Mh-  At  first  glance  it  seems  as 
though  nothing  has  been  gained  by  including  the  halo  component.  The  rotation  curve 
of  the  Plummer  sphere  declines  and  approaches  the  Keplerian  curve  also  at  large  radii. 
There  is  a crucial  difference,  however,  between  the  halo  and  the  disk.  Since  the  disk  is 
luminous,  its  exponential  scale  length  can  be  observed.  The  halo,  though,  is  made  up 
of  dark  matter  which  is  not  accessible  to  direct  observation.  The  totality  of  knowledge 
about  the  matter  distribution  in  the  halo  must  be  inferred  by  its  gravitational  effects. 
The  end  result  is  that  the  parameter  b,  unlike  the  exponential  disk  scale  length  eo,  is 
free,  constrained  only  by  the  observed  rotation  curve.  Hence  the  value  of  b can  be  near 
or  beyond  the  observed  disk  (optical  or  HI),  thereby  placing  the  disk  completely  within 
the  rising  part  of  the  halo  rotation  curve. 

For  modeling  the  bar,  the  triaxial  homeoidal  density  distribution  given  by 


™<;  (2-43) 

to  772  > 1, 

O 7/^  2 2 

where  = fr  + fr  + TT  and  a > & > c > 0 are  the  long,  intermediate,  and  short  bar 

Ct«  o 

axes,  respectively,  is  used.  This  is  the  n = 2 Ferrers  potential  (Ferrers  1877;  see  also 
Binney  and  Tremaine  1987,  p.61).  This  particular  bar  component  was  chosen  not  only 
because  it  represents  many  features  of  observed  bars  rather  well — Hunter  et  al.  (1988) 
and  Ball  (1984,  1992)  found  that  bar  components  with  a Gaussian  brightness  profile, 
similar  to  that  of  the  n - 2 Ferrers  bar,  best  fit  the  infrared  observations  of  the  bars 
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of  NGC  3992  and  NGC  3359,  respectively — but  also  because  its  dynamics  has  been 
the  subject  of  a number  of  previous  investigations  (e.g.  de  Vaucouleurs  and  Freeman 
1972;  Papayannopoulos  and  Petrou  1983;  Athanassoula  et  al.  1983;  Pfenniger  1984a, 
1985).  The  strict  inequalities  in  the  bar  axes’  lengths  are  only  required  for  mathematical 
reasons  in  the  derivation  of  the  associated  potential  and  forces.  They  do  not  represent 
a real  restriction  from  essentially  spherical,  oblate  or  prolate  figures.  Exact  solutions 
for  these  special  cases  already  exist. 


Figure  2^:  Rotation  curves  of  the  Plummer  sphere  (solid 
curve)  and  a point  with  the  same  total  mass  (dotted  curve). 

The  central  density  pc  of  the  bar  can  be  expressed  in  terms  of  the  total  mass  Mp 

and  its  axes’  lengths  as 


105  Mb 


pc 


32tt  ahc 


(2-44) 
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The  projected  surface  density  of  this  figure  is  obtained  by  integrating  p(x,y,z)  with 


respect  to  z over  the  range  where  p is  nonzero: 

Zrnax 


p{x,y,z)dz, 


— Z 


max 


where  Zmax  — / 1 


jf.  The  result  is  expressed  by 


^B{x,y)  = ^cpc 
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(2-45) 


(2-46) 


The  potential  in  the  disk  plane  can  be  written  as 

VB{x,y)  = -^mWnoxV  - WnoxS/  - W2wxV 

D 


kkiooy 


+lTo2oa:^  + 1^2002/^  — kkoio^:^) 


(2-47) 


+Vfooo  — — VV3002/^}? 

where  C = 27rGpcabc  = ^GMb,  and  the  Vkp  are  coefficients  dependent  on  both  the 
axes’  lengths  a,  b,  and  c and  the  position  {x,y).  Complete  derivations  of  the  general 
three-dimensional  form  of  the  potential,  the  Cartesian  components  of  the  force,  and  the 
coefficients  ITj/yt  are  given  by  Pfenniger  (1984a,  appendix).  A summary  of  all  the  bar 
quantities  needed  for  this  study  is  included  in  the  appendix  of  this  thesis. 

Despite  the  fact  that  there  is  no  “axisymmetric”  rotation  velocity  associated  with 
this  bar  potential,  an  estimate  of  this  quantity  can  be  obtained.  By  analogy  with  the 
axisymmetric  disk  and  halo  rotation  velocities,  the  “circular”  velocity  VB{r,0)  due 
solely  to  the  bar  potential  is  given  by 


VB(r,0) 


JVB{r,6) 

dr 


(2-48) 


Above,  the  velocity  and  potential  have  been  taken  in  terms  of  the  polar  coordinates  (r,9) 
instead  of  the  Cartesian  ones  (x,y).  The  “axisymmetric”  rotation  velocity  u^(r),  then,  is 

7T 

2 / 


VB{r)  = 


2 

7T 


or 


(2-49) 
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In  a similar  manner,  all  other  “axisymmetric”  quantities  of  the  bar  can  be  estimated. 
Figure  2-5  gives  representative  “axisymmetric”  rotation  curves  for  three  sets  of  bar 
dimensions. 


Figure  2-5;  Three  nomralized  “axisymmetric”  bar  rotation  curves  for  cases  where  a = 

1 and  c = 0.1:  (1)  b = 0.15,  (2)  b = 0.3,  and  (3)  b = 0.6.  The  dotted 
line  represents  the  Keplerian  rotation  curve  for  a point  mass  equal  to  Mg. 

As  in  the  original  CG  method,  the  observed  angle-averaged  rotation  curve  is  fit 

by  the  model  rotation  curve.  The  model  curve  in  this  case,  though,  is  derived  from  the 

component  rotation  curve  by  the  relation 

v'^{r)  = v]){r)■^rv]J{r)  + v%{r),  (2-50) 

where  vj'{r)  is  the  total  rotation  velocity.  The  method  of  curve  fitting  is  somewhat  more 
crude  in  our  modified  method  than  in  the  original  CG  method.  In  particular,  the  use 
of  a least  squares  fitting  routine  is  forgone  in  lieu  of  simple  visual  fitting.  There  are  a 
couple  of  reasons  for  this.  First,  in  barred  spiral  galaxies  the  noncircular  motions  can  be 
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much  more  pronounced  than  in  normal  spirals  (e.g.  Duval  and  Athanassoula  1983).  The 
uncertainty  in  the  “correct”  axisymmetric  rotation  curve,  then,  grows  in  proportion  to 
the  magnitude  of  the  noncircular  motions.  Large  uncertainties,  in  themselves,  however, 
do  not  preclude  the  use  of  a mathematical  fitting  scheme.  A second  and  more  serious 
problem  was  that  least  square  fits  in  which  all  input  parameters  were  allowed  to  adjust 
freely  resulted,  in  some  cases,  in  physical  implausibilities,  such  as  near- zero  bar  masses. 
For  these  reasons  the  following  procedure  for  fitting  rotation  curves  was  implemented. 

First,  since  the  component  contributions  to  the  rotation  curve  are  somewhat 
independent  of  each  other  (the  bar  component  dominates  the  inner  region  of  the  observed 
rotation  curve,  the  disk  component  the  central  region,  and  the  halo  component  the  outer 
region),  it  was  decided  to  fit  visually  one  component  at  a time,  starting  from  the  center 
outward.  The  first  major  contributor  to  the  rotation  curve  is  the  bar.  The  bar  that  was 
used  was  the  one  with  the  maximum  allowable  mass  under  the  constraints  of  the  rotation 
curve.  This  component  typically  fit  the  steep  initial  rise  of  the  rotation  curve.  The  next 
component  included  was  the  exponential  disk.  Again,  the  maximum  disk  allowable  by 
the  rotation  curve  was  used.  The  addition  of  the  disk  component  permitted  good  fits  to 
the  observed  rotation  curve  over  most  of  the  radial  range.  The  Plummer  sphere  halo 
component  was  added  last  in  order  to  “hold  up”  the  rotation  curves  at  large  radii.  The 
parameter  values  obtained  this  way  were  adjusted  slightly  in  order  to  achieve  the  best 
visual  fit.  The  parameters  that  were  adjusted  in  order  to  provide  the  fits  to  the  rotation 
curve  are  given  in  Table  2-3.  While  the  bar  contribution  to  the  rotation  curve  certainly 
depends  on  the  axes’  lengths  a,  b,  and  c,  these  parameters,  as  we  will  see,  are  fixed 


by  the  surface  photometry. 
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Table  2-3:  Adjustable  parameters  in  the  axisymmetric  rotation  curve  fitting  procedure. 

Component  Adjustable  Parameter(s) 

Ferrers  Bar  Mb 

Exponential  Disk  cq  and  eo 

Plummer  Sphere  Halo  M//  and  b 

The  method  of  Stark  (1977)  is  used  to  derive  the  lengths  of  the  bar’s  axes.  A 
complete  explanation  of  the  method  can  be  found  in  Stark  (1977)  and  in  Ball  (1984), 
thus  the  derivation  will  not  be  included  here.  In  short,  Stark’s  method  assumes  that  the 
volume  brightness  of  the  bar  is  constant  on  similar,  nested,  triaxial  ellipsoidal  surfaces. 
Contopoulos  (1956)  has  shown  that  the  projected  isophotes  of  such  a bar,  from  any 
viewing  angle,  are  similar,  concentric  ellipses.  Assuming  that  the  triaxial  bar  figure 
has  one  of  its  principle  axes  normal  to  the  disk  plane,  one  can  combine  the  observed 
axial  ratio  of  the  bar  isophotes,  /?o,  the  inclination  angle  of  the  galaxy,  i,  and  the  angle 
between  the  line  of  nodes  of  the  disk  and  the  major  axis  of  the  bar  isophotes,  tp,  to 
give  a one  parameter  family  of  triaxial  figures  as  the  solution.  The  parameter  <f>,  which 
controls  the  family,  is  in  some  sense  a measure  of  the  true  azimuth  of  the  triaxial  figure 
in  the  disk  plane  referred  to  some  fiducial  direction.  In  the  present  analysis  p is  taken 
to  be  the  angle  measured  in  the  disk  plane  between  the  long  axis  of  the  bar  and  the 
line  of  nodes  of  the  disk.  A useful  property  of  this  family  of  solutions  is  that  there 
is  a one-to-one  correspondence  between  p and  the  axial  ratios  of  the  bar.  Specifying, 
say,  the  axial  ratio  da  is  equivalent  to  specifying  (f>  directly,  and  that  is  the  actual 
procedure  used  here.  The  value  of  da  is  taken  to  be  1/10  after  Pfenniger  (1984a  [see 
also  references  therein,  de  Vaucouleurs  and  Freeman  1972;  Kormendy  1982],  1985).  In 
all  of  the  cases  considered  in  this  thesis  the  long  and  intermediate  bar  axes  a and  b are 


assumed  to  lie  in  the  principle  disk  plane  of  the  galaxy,  and  the  short  axis  c is  assumed 
to  be  perpendicular  to  the  disk  plane. 


After  the  determination  of  the  lengths  of  the  bar  axes  from  the  Stark  analysis,  and 
the  visual  fitting  of  the  axisymmetric  rotation  curve,  we  derive  starting  values  for  Qp 
and  for  the  parameters  which  define  the  spiral  perturbation.  There  is  ample  theoretical 
and  numerical  evidence  (Contopoulos  1981;  Schwarz  1981;  Teuben  and  Sanders  1985; 
Sellwood  1980,  1981;  Thielheim  and  Wolff  1981,  1984;  Sparke  and  Sellwood  1987) 
that  bars  end  at  or  slightly  inside  corotation,  the  distance  at  which  the  angular  rotation 
rate  of  matter  in  the  galaxy  equals  flp.  For  a first  estimate  of  the  pattern  speed,  then, 
the  bar  is  assumed  to  end  at  corotation.  As  before,  the  parameters  characterizing  the 
spiral  pattern  beyond  the  bar  are  derived  from  the  available  surface  photometry.  A 
new  feature  is  now  added,  however,  to  the  form  of  the  perturbation.  Since  the  spiral  is 
assumed  to  begin  at  or  near  the  end  of  the  bar,  an  inner  cutoff  is  added  to  the  formula 
for  the  spiral  amplitude  Amir): 


As  we  shall  see  later,  the  optimum  positions  for  these  inner  cutoffs  yield  important 
information  concerning  the  nature  of  the  bars  in  the  most  successful  models.  Also, 
some  changes  to  the  argument  of  the  cosine  in  Eq.  (2-6)  are  required.  Namely, 


(2-51) 


argument  = < 

At  this  point,  the  entire  barred  spiral  model  is  specified. 


r > a 


(2-52) 


r < a. 


As  in  the  original  CG  method,  the  next  step  is  to  calculate  the  main  families 


of  periodic  orbits.  In  the  barred  spiral  models,  periodic  orbit  families  were  searched 
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for  between  the  ILR  and  the  8/1  resonance  inside  corotation,  and  between  the  -8/1 
and  -2/1  (outer  Lindblad)  resonance  outside  corotation.  For  bars  of  moderate  to  high 
strength,  however,  there  are  substantial  regions  of  “energy”  r^,  both  inside  and  outside 
of  corotation,  where  there  are  no  stable,  dynamically  significant  periodic  families.  In 
order  to  represent  matter  existing  in  these  regions  for  the  surface  response  density 
calculation,  we  included  orbits  started  here  on  “circular”  trajectories.  The  initial  radius 
of  a “circular”  orbit  of  energy  rc  was  taken  to  be  r^,  and  the  initial  radial  velocity 
was  taken  to  be  zero  (hence  the  designation  “circular”).  Since  the  majority  of  orbits  in 
these  regions  turn  out  to  be  stochastic,  the  surface  density  response  they  generate  is  not 
particularly  sensitive  to  the  particular  choice  of  initial  conditions.  Thus  we  have  given 
them  the  most  simple  initial  conditions  possible,  the  “circular”  orbit  initial  conditions 
described  above.  Despite  the  fact  that  most  of  these  orbits  turn  out  to  be  stochastic, 
they  yield  interesting  density  responses  nonetheless,  as  we  will  see  later  in  Chapters 


3,  4,  and  5. 


Next,  the  surface  response  density  map  is  created.  The  procedure  is  exactly  the 
same  as  in  the  original  CG  method,  except  for  the  handling  of  the  “circular”  orbits 
described  above.  The  differences  are  (1)  in  the  formula  for  the  weight  Wm,  which 
gauges  the  relative  amount  of  trapping  done  by  the  m*  family,  Xm,  the  deviation  of  the 
orbit  from  r^,  is  taken  to  be  1 kpc,  a typical  value  for  this  type  of  orbit,  (2)  the  starting 
points  for  dispersed  orbits  along  the  “circular”  orbits  are  taken  to  be  every  interval 
2xNOA  azimuth,  (3)  in  cases  where  the  given  “energy”  and  initial  conditions  of  a 
“circular”  orbit  result  in  an  imaginary  angular  momentum  solution  (i.e.  where  the  initial 
conditions  lie  outside  the  limiting  curve)  the  “energy”  is  artificially  increased  to  the  point 
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where  the  angular  momentum  becomes  real  (this  procedure  is  justified  by  noting  that 
both  observations  of  real  barred  galaxies  and  N-body  simulations  of  barred  galaxies 
show  that  these  regions  are  populated  with  stars),  and  (4)  all  “circular”  orbits  and  orbits 
dispersed  around  them  are  integrated  for  a fixed  time  T,  chosen  such  that  good  estimates 
of  the  average  density  response  are  derived  from  these  orbits  and  computation  times  axe 
held  to  reasonable  lengths.  Typically  T is  chosen  to  be  approximately  1 billion  years. 

These  changes  require  some  changes  in  the  precise  form  of  the  weighting  factors; 
however,  the  meaning  of  each  factor  remains  the  same.  Specifically,  for  all  orbits, 

77 

2 

ms  = Co  exp  (-eoxc)  + “ / ^^(rc,  6)dd,  (2-53) 

0 

and  for  the  “circular”  orbits, 

1 

^time  — (2—54) 

The  resulting  density  map  is  preprocessed  exactly  as  before  prior  to  the  analysis  for 
self-consistency. 

The  essence  of  the  analysis  of  the  density  response  map  remains  unchanged  in  the 
modified  method  (i.e.  the  application  of  the  CG  method  to  barred  spiral  galaxies); 
however,  one  detail  concerning  the  computation  of  the  amplitudes  of  the  Fourier 
components  of  the  imposed  density  has  been  altered.  Since  the  bar  potential  is  not 
put  into  the  model  explicitly  in  the  form  of  its  Fourier  components,  the  amplitude  of 
the  29  component  of  the  imposed  density  is  derived  in  exactly  the  same  way  as  the 
29  component  of  the  response  density;  specifically,  the  imposed  density  is  calculated 
as  a function  of  azimuth,  for  a fixed  semiannulus  radius,  at  exactly  the  same  azimuthal 
resolution  as  the  response  density  map.  The  imposed  density  at  a given  position  is 
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derived  from  (1)  the  surface  density  of  the  exponential  disk  and  the  projected  bar, 
and  (2)  the  spiral  perturbation  according  to  the  formula  of  Vandervoort.  A Fast  Fourier 
Transform  is  then  performed  on  this  array  of  imposed  densities  at  the  same  time  as  on  the 
array  of  response  densities.  At  each  radius,  then,  the  amplitudes  of  the  20  components 
are  compared  by  the  ratio  of  response  to  imposed,  as  in  the  original  method.  The  phase 
difference  calculation  is  unchanged. 

Gas  Response  Using  Smoothed  Particle  Hydrodynamics  (SPH) 

To  complement  the  self-consistent  stellar  models  achieved  by  the  modified  method, 
the  gas  response  of  the  “best”  model  in  each  case  has  been  investigated  via  the  use  of 
a two-dimensional  smoothed  particle  hydrodynamics  (SPH)  code  generously  provided 
to  the  author  by  Dr.  Nikos  Hiotelis.  The  SPH  method  was  introduced  independently 
by  Lucy  (1977)  and  Gingold  and  Monaghan  (1977)  and  uses  a Monte  Carlo  method 
to  solve  the  equations  of  hydrodynamics.  The  complete  details  of  the  SPH  method  are 
given  in  several  review  articles  (Hemquist  and  Katz  1989;  Benz  1990;  Monaghan  1985; 
Steinmetz  and  Muller  1992)  and  references  therein,  thus  only  a brief  overview  of  the 
general  method  (in  two  dimensions)  will  be  given  here,  along  with  the  relevant  details 
of  the  specific  SPH  implementation  used  for  the  present  calculations. 

SPH  is  a Lagrangian  method  which  does  not  require  the  use  of  a computational 
grid;  therefore,  unlike  finite-difference  schemes,  which  require  a grid,  no  computational 
resources  are  wasted  simulating  large  voids  that  may  arise  (Steinmetz  and  Muller  1992). 
Since  the  SPH  method  is  not  very  complicated,  it  does  not  require  a lengthy  code; 
furthermore,  it  is  quite  robust.  Steinmetz  and  Muller  (1992)  mention  also  three  new 
developments  which  have  greatly  enhanced  the  potential  of  SPH:  (1)  the  achievable 
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resolution  has  been  improved  dramatically  with  the  advent  of  the  variable  smoothing 
length,  (2)  global  time  step  restrictions  have  been  alleviated  via  the  use  of  a separate  time 
step  for  each  particle,  and  (3)  the  introduction  of  the  hierarchical  tree  method  to  calculate 
the  gravitational  forces  has  reduced  computation  times  for  SPH  with  self-gravity  from 
the  order  of  to  MogA^  without  altering  its  original  Lagrangian  formulation.  Finally, 
the  SPH  scheme  and  hierarchical  tree  method  can  be  vectorized  to  take  advantage  of 
highly  parallelized  computing  architectures. 

In  the  SPH  method  a continuous  fluid  medium  is  modeled  as  an  ensemble  of  N 
fluid  elements  (in  this  thesis  N is  approximately  12,000  in  all  cases),  each  taken  to  be 
a “particle”  with  position  r,-  and  mass  mi,  smoothed  out  according  to  some  chosen 
smoothing  kernel  W{r  - r'.  A).  The  kernel  is  a function  strongly  peaked  around 


r — r 


= 0,  and  the  quantity  A is  called  the  smoothing  length.  The  kernel  must 


fulfull  certain  requirements.  First,  it  must  guarantee  conservation  of  momentum,  both 
linear  and  angular.  An  easy  way  to  do  this  is  to  make  the  kernel  spherically  symmetric. 


Secondly,  it  must  have  the  properties 


J W{r,\)d‘^r  = 1 


and 


lim  mr,A)  = (^(r). 
A-^O 


(2-55) 


The  ensemble  of  particles  together  with  the  smoothing  kernel  can  be  used  to  estimate 
the  mean  value  of  any  spatially  varying  physical  variable  at  any  arbitrary  position  r. 


The  actual  mean  value  at  r of  the  spatially  varying  function  A(r)  is  given  by 


(A(r))  = J A(r')W (r  — r' , X)(Pr' . 


(2-56) 


The  estimate  using  the  “smoothed”  particles  is  given  by  the  sum 


N 


(^(r))  = ^mj-^lT"(r- rj,A), 

PAj) 


j=i 


(2-57) 
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where,  for  the  density  p(Xj),  the  similarly  estimated  value 

N 

2 = 1 


(2-58) 


is  used.  The  time  derivative  and  gradient,  respectively,  of  the  quantity  (A(r))  can  be 


written  (Steinmetz  and  Muller  1992)  as 

d 


and 


(2-59) 


V(A(r))  = j A{v')\7W{r-v',X)d^r'. 


These  results  can  be  used  to  derive  the  averaged  hydrodynamic  equations  of 

motion.  Specifically,  they  are  obtained  by  multiplying  each  term  of  the  exact  equations 

by  the  smoothing  kernel  and  summing  over  all  particles  (e.g.  Monaghan  and  Gingold 

1983;  Monaghan  1985;  Benz  1990).  Including  also  the  artificial  viscosity  terms,  the 

following  set  of  equations  (Hiotelis  et  al.  1991)  are  obtained: 

dvi 


dt 


dui 

dt 


1 

2 


N 

j=i 

>«•  , p/ 
[p1  p)_ 

(1  + Uij)ViW 

dvi  _ 

V,', 

dt 

f 7 

'P^  , Pi 
/i  !>)  J 

(1  + n,y)(vi  - 

(2-60) 


Above,  Vf,  Pi,  pi,  Ui,  Vi,  and  Vj-l-T  are  the  velocity,  pressure,  density,  thermal  energy 
per  unit  mass,  gravitational  potential,  and  gradient  of  the  kernel  W,  respectively,  at  the 
position  T;  of  the  particle.  The  term  II, y includes  the  effect  of  artificial  viscosity. 
The  form  of  II, y is  taken  from  Monaghan  and  Gingold  (1983)  and  is  the  same  as  that 
used  by  Hiotelis  et  al.  (1991): 


n 


r\ 

Oipij  + PPij- 


(2-61) 
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Here,  a and  (3  are  constant  coefficients  and  j.iij  is  defined  by 

l^ij  = i cij(k.-rfi  W) 

0,  otherwise. 


(2-62) 


Above,  Ay  and  cij  are  the  average  smoothing  length  and  sound  speed,  respectively. 


y 


between  positions  r,-  and  r^'.  The  parameter  = 0.01  A^  is  introduced  to  avoid  possible 
numerical  divergences  in  mj.  In  the  expression  for  Ily,  the  term  linear  in  /ty  produces 
a shear  and  bulk  viscosity  (Monaghan  1985),  while  the  quadratic  term  is  included  to 
handle  high  Mach  number  shocks  and  is  roughly  equivalent  to  the  Von  Neumann- 
Richtmyer  viscosity  used  in  finite-difference  methods.  In  all  calculations  done  here,  a 
and  /?  are  taken  to  be  unity.  According  to  Hernquist  and  Katz  (1989),  these  values 
represent  a fair  compromise  for  the  typical  range  of  Mach  numbers  encountered. 

In  order  to  integrate  the  hydrodynamic  equations  an  equation  of  state  must  be 
specified.  The  gas  in  this  case  is  assumed  to  behave  like  an  ideal  gas.  Therefore,  the 
pressure,  density  and  thermal  energy  of  the  gas  is  related  by 


Pj  = (7  - l)piUi,  (2-63) 

where  7 is  defined  to  be  the  ratio  of  specific  heats.  Since  the  calculations  are  done 
in  two  dimensions  instead  of  three,  the  values  of  7 corresponding  to  the  adiabatic  and 
isothermal  cases  are  2 and  1 , respectively.  The  actual  behavior  of  interstellar  gas  most 
probably  is  somewhere  between  these  two  extremes,  therefore  a compromise  value  of  7 
= 5/3  is  used.  The  wide  variation  present  in  the  temperature  of  the  interstellar  medium 
presents  a problem  concerning  the  appropriate  choice  of  the  thermal  energy  ui.  An 
initial  value  of  w/  corresponding  to  an  average  temperature  of  around  3500  K is  used. 
We  consider  this  a fair  compromise  between  the  larger  expanses  of  very  hot,  low  density 
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interstellar  gas  (T  ft;  10®  K,  n ft;  10“^  cm~^)  and  the  smaller  expanses  of  relatively 
cooler  and  denser  HI  clouds  (10^  K < T < 6 x 10^  K,  20  cm“^  > 0.3  cm~^)  and 
HII  regions  (T  ft^  8 x 10^  K,  n > 0.5  cm”^). 

It  still  remains  to  specify  the  exact  smoothing  kernel  used.  The  two-dimensional 
version  of  the  spherically  symmetric  three-dimensional  spline-based  kernel  proposed  by 
Monaghan  and  Lattanzio  (1985)  and  employed  by  Hiotelis  et  al.  (1991)  is  used.  It 
is  defined  by 


1 - -h  0 < u < 1 

i(2-^^)^  l<u<2 

0,  u > 2, 


(2-64) 


where  u = I . Also,  the  smoothing  length  is  assumed  to  be  variable  in  space  and  time, 
hence  A = A(r,  t).  This  is  done  to  increase  the  resolution  possible  with  the  SPH  method, 
since  A directly  limits  the  resolution  (all  structure  on  scales  smaller  than  A are  strongly 
smoothed  out).  Whereas  the  introduction  of  a time-dependent  A is  straightforward,  there 
exist  two  possibilities  to  introduce  a spatially  variable  A (Steinmetz  and  Muller  1992). 
These  are  the  “gather”  and  “scatter”  approaches  (cf.  Hernquist  and  Katz  1989).  In  the 
“gather”  method,  the  smoothed  average  value  of  the  arbitrary  physical  quantity  A(r) 
is  calculated  by  summing  the  contributions  from  the  whole  of  space  weighted  by  the 
smoothing  kernel  centered  on  r using  the  local  value  A(r)  of  the  smoothing  length: 


{A{r))  gather  = / A(r')  fP  (r  - r,  A(r))  dV' . (2-65) 


Hence  the  contributions  to  A at  r are  “gathered”  in  from  the  surrounding  elements, 
weighted  appropriately  by  the  kernel  at  r.  The  “scatter”  method,  on  the  other  hand, 
derives  its  estimate  for  the  same  quantity  by  using  the  kernel  local  to  the  contributing 
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element  and  “scattering”  its  appropriate  contribution  to  position  r: 


{A{r)) scatter  = A{v')W {v  - V , \{r'))(fr' . (2-66) 


Exclusive  use  of  one  or  the  other  of  these  approaches  does  not  guarantee  the  conservation 
of  energy,  linear  momentum,  or  angular  momentum.  A combination  is  usually  used. 
Benz  (1990)  has  symmetrized  the  smoothing  length  by  using  Wij  = W(ri  - rj,  (Aj  + Xj)  / 
2),  while  Hernquist  and  Katz  (1989)  have  symmetrized  the  kernel  itself  Wij  = 0.5(Vk(r/ 
- rj.  A;)  + W{Vi  - Tj,  Xj)).  In  the  present  work,  the  first  method  of  symmetrizing  the 
smoothing  length  is  used.  Two  further  questions  remain  concerning  A(r,  t):  (1)  How 
should  A be  defined?  and  (2)  how  should  it  evolve  during  the  integration? 

It  seems  natural  to  associate  the  smoothing  length  with  the  local  density  such  that 
denser  regions  have  smaller  smoothing  lengths  in  order  to  provide  better  resolution. 
One  way  to  do  this  (Steinmetz  and  Muller  1992)  is  to  let 

Vr.)  = c(^)',  (2-67) 

where  C is  a parameter  of  order  unity.  The  problem  with  this  definition  is  that  A(r,)  is 
needed  in  order  to  compute  p(ri)  in  the  first  place.  One  way  around  this  problem  is  to 
use  the  value  of  the  density  at  the  previous  time  step  to  calculate  the  smoothing  length 
(Miyama  et  al.  1984).  Another  way,  proposed  by  Benz  (1990),  is  to  take  the  time 
derivative  of  both  sides  above  and  eliminate  the  density  altogether  using  the  equation 
of  continuity  to  yield 


dX[v) 

dt 


1 

-AV  • V. 
2 


(2-68) 


In  this  way  the  smoothing  length  can  be  considered  to  be  just  another  hydrodynamic 
variable  to  be  evolved  in  the  course  of  the  integration.  Yet  another  method  (Hernquist 
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and  Katz  1989)  is  to  link  the  value  of  A to  the  number  of  neighboring  particles  (within 
a radius  of  2A),  and  let  A change  in  space  and  time  so  as  to  keep  this  number  fixed.  The 
method  used  here  is  the  second  one  where  A is  given  an  initial  value  for  each  particle 
and  is  subsequently  evolved  according  to  Eq.  (2-68)  above. 

The  algorithm  used  to  integrate  the  hydrodynamic  equations  of  motion  is  a 
predictor-corrector  time  stepping  rule.  At  any  given  point  in  the  integration  a single 
time  step  St  is  used  to  advance  all  the  “particles.”  The  stepsize  varies  from  one  step  to 
another  and  is  given  by  the  relation 


St  — min 


(2-69) 


is  the  estimate  of  the  circular  orbit  period  at  position  r,-,  and 
(the  Courant  number)  is  an  input  parameter  of  order  0.5  utilized  to  suppress  numerical 
instabilities.  In  all  cases  we  take  c/v  to  be  0.4.  Also,  F,-  and  v,-  are  the  acceleration  and 
velocity  of  the  particle.  The  terms  in  the  expression  for  St  are  included  to  allow, 
respectively,  fast-moving  “particles,”  highly  accelerated  “particles,”  and  “particles”  very 
close  to  the  center  to  be  followed  accurately. 


where  T — 2Tr 


CHAPTER  3 
PROGRAM  GALAXIES 


The  galaxies  chosen  for  this  study  are  NGC  3992,  NGC  1073,  and  NGC  1398, 
and  are  shown  in  Figures  3-1,  3-2,  and  3-3.  These  galaxies  were  chosen  for  several 
reasons.  First,  the  neutral  hydrogen  emission  from  all  three  galaxies  has  been  observed 
at  the  Very  Large  Array  (VLA)  of  the  National  Radio  Astronomy  Observatory  (NRAO) 
at  high  resolutions  and  with  high  signal-to-noise  ratios.  These  observations  reveal 
information  about  the  kinematics,  and  therefore  the  total  gravitational  potential,  of  the 
disk.  Secondly,  detailed  surface  photometry  exists  for  each  of  the  galaxies.  Elmegreen 
and  Elmegreen  (1985)  have  provided  both  blue  and  near-infrared  surface  photometry  of 
NGC  3992  and  NGC  1073.  Ohta  et  al.  (1990)  have  observed  NGC  1398  in  the  B band. 
In  addition,  Grosbpl  (1985)  has  scanned  and  analyzed  the  images  of  all  three  galaxies 
on  the  red  Palomar  Sky  Survey  copy  plates.  The  surface  photometry  yields  information 
concerning  the  distribution  of  luminous  matter  in  the  disk.  Finally,  these  galaxies  span 
a moderately  wide  range  of  barred  spiral  morphological  types,  de  Vaucouleurs  et  al. 
(1976)  classify  NGC  1398  as  (R')SB(r)ab,  NGC  3992  as  SB(rs)bc,  and  NGC  1073  as 
SB(rs)c.  This  latitude  allows  us  to  look  for  possible  correlations  between  the  model 
parameters  of  the  most  successful  models  and  morphological  type.  Furthermore,  since 
these  galaxies  show  no  gross  peculiarities,  they  can  be  considered  to  be  representative 
examples  of  their  respective  types. 
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NGC  3992 


Observations 

Neutral  Hydrogen 


High-resolution  neutral  hydrogen  observations  of  NGC  3992  were  obtained  by 
Gottesman,  Ball,  Hunter,  and  Huntley  (1984,  hereafter  GBHH)  at  the  VLA  between 
September,  1980  and  January,  1983.  The  effective  resolution  obtained  was  26."  1 x 20. "0, 
with  the  major  axis  of  the  Gaussian  beam  at  a position  angle  of-22.°l.  In  comparison, 
the  diameter  of  the  HI  disk  is  504"  and  the  position  angle  of  the  line  of  nodes  is  -1 1 1.°5 
± 0.°6.  The  full- width  half-power  (FWHP)  velocity  resolution  was  25.2  km  s~^  The 
reader  is  referred  to  GBHH  for  further  details  of  the  data  reduction. 

The  angle-averaged  rotation  curve  was  derived,  along  with  the  following  orienta- 
tion parameters  for  the  HI  disk,  via  a least-squares  fit  to  the  velocity  field:  major  axis  at 
a position  angle  of  -lll.°5  + 0.°6,  inclination  of  53.°4  + 0.°9,  and  a systemic  velocity 
of  1045.8  + 0.6  km  s“^.  The  data  points  so  derived  are  given  in  Table  3-1.  The  quoted 
errors  are  the  formal  errors  of  the  least-squares  solution.  The  angular  radii  of  columns 
1 and  4 in  Table  3-1  have  been  converted  to  linear  dimensions  in  columns  2 and  5 using 
an  assumed  distance  to  NGC  3992  of  14.2  Mpc  (de  Vaucouleurs  1979).  The  data  in 
Table  3-1  are  plotted  as  filled  circles  with  error  bars  in  Figure  3-4. 

GBHH  find  that  the  spiral  structure  of  NGC  3992  is  only  weakly  discernible  in 
the  neutral  hydrogen  distribution.  The  optical  image,  nevertheless  (Figure  3-2),  is  seen 
clearly  to  reveal  spiral  structure  which  extends  from  an  incomplete  ring  around  the  bar. 
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Table  3-1:  HI  rotation  velocity  data  for  NGC  3992. 


R (arcmin) 
(1) 

R (kpc) 
(1) 

V (km  s 
(2) 

R (arcmin) 
(1) 

R (kpc) 
(1) 

V (km  s ^) 
(2) 

0.14 

0.58 

141.8  + 37.8 

2.89 

11.9 

269.1  + 5.5 

0.37 

1.53 

123.6  ± 52.7 

3.13 

12.9 

270.9  + 5.5 

0.64 

2.64 

170.9  ± 32.7 

3.36 

13.9 

267.3  ± 7.1 

0.84 

3.47 

195.6  ± 16.5 

3.63 

15.0 

265.5  ± 7.5 

1.14 

4.71 

222.5  ± 10.9 

3.86 

15.9 

258.2  + 10.9 

1.37 

5.66 

240.0  ± 9.1 

4.09 

16.9 

249.1  ± 11.1 

1.64 

6.77 

254.2  + 7.3 

4.36 

18.0 

247.3  ± 12.5 

1.90 

7.85 

258.2  + 5.8 

4.63 

19.1 

240.0  ± 10.9 

2.14 

8.84 

262.2  ± 4.7 

4.89 

20.2 

240.0  + 20.0 

2.36 

9.75 

265.5  ± 4.2 

5.10 

21.1 

238.2  ± 18.2 

2.63 

10.9 

267.3  + 4.4 

(1)  Radius  measured  from  the  center  of  the  galaxy 

(2)  Rotation  velocity 

They  also  find  a deficiency  of  HI  inside  this  optical  ring.  Two  possibilities  exist.  Either 
the  total  gas  surface  density  suffers  a real  depression  in  the  central  region  of  the  galaxy, 
or  most  of  the  hydrogen  in  this  region  has  been  converted  to  molecular  form  and  is 
not  being  seen  in  the  HI  observations.  GBHH  cite  observations  of  NGC  3992 
by  Young,  which  set  an  upper  limit  to  the  H2  surface  density  of  only  a few  times 
larger  than  the  neutral  hydrogen  surface  density,  and  conclude  that  the  total  gas  surface 
density  near  the  center  appears  to  be  less  than  in  the  disk.  They  add  that  “it  is  possible 
that  dynamical  effects  associated  with  the  bar  may  have  exaggerated  the  phenomenon” 
(GBHH,  pp.  477-478).  As  we  shall  see  from  the  results  of  our  gas  calculations  (Chapter 
6),  bars  can  indeed  “sweep  out”  gas  from  large  regions  inside  the  bar  radius. 
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Figure  3-4:  Comparison  of  observed  and  theoretical  rotation  curves  for  NGC  3992.  The 
theoretical  curve  is  derived  from  our  most  successful  model  of  NGC  3992.  Also  shown  are  the 
contributions  of  the  separate  components  of  this  model  to  the  total  theoretical  rotation  curve. 

As  far  as  the  kinematics  of  NGC  3992  are  concerned,  GBHH  find  that  the  motion 
of  the  gas  is  dominated  by  circular  rotation,  with  only  weak  irregularities  associated 
with  the  optical  spiral  arms.  They  interpret  these  irregularities  as  streaming  of  the  gas 
along  the  spiral  arms.  Also  noted  is  a discontinuity  in  the  isovelocity  contours  of  the 
velocity  field  across  the  major  axis  of  the  galaxy,  which  the  authors  tentatively  interpret 
as  an  effect  of  the  flow  of  gas  around  the  bar  of  NGC  3992.  The  uncertainty  in  this 
interpretation  stems  from  the  limited  resolution  in  this  region  due  to  low  gas  densities. 
Finally,  GBHH  note  some  peculiarities  in  the  neutral  hydrogen  velocities  along  the 
major  axis  but  at  radial  distances  greater  than  about  3. '8  (15.7  kpc).  Specifically,  the 
velocities  show  a rather  sudden  decrease  at  this  radius,  followed  by  a more  gradual 
decline  thereafter.  GBHH  identify  this  kinematical  feature  with  a possible  truncation 


of  the  disk. 
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Table  3-2:  Selected  global  and  disk  parameter  values  adopted  for  NGC  3992.  The  errors 
given  are  simply  the  formal  errors  of  a least-squares  analysis  of  the  observed  velocity  field 
and  do  not  imply  that  these  quantities  have  been  determined  to  this  level  of  precision. 


Parameter 

Value 

Systemic  velocity  (heliocentric,  km  s“^)^ 

1045.8  + 0.6 

Distance  (Mpc)'’ 

14.2 

Inclination  angle  (°)® 

53.4  + 0.9 

Position  angle,  line  of  nodes  (°)® 

-111.5  + 0.6 

^Gottesman,  Ball,  Hunter,  and  Huntley  (1984) 
'’de  Vaucouleurs  (1979) 

Surface  Photometry 

Elmegreen  and  Elmegreen  (1985,  hereafter  EE)  have  provided  blue  {B  passband) 
and  near-infrared  (/  passband)  surface  photometry  of  fifteen  barred  spiral  galaxies, 
including  NGC  3992.  They  photographed  these  galaxies  between  1979  and  1981  using 
the  1.2  meter  Palomar  Schmidt  telescope.  The  blue  images  were  taken  using  baked 
103a-O  emulsions  with  a GG  385  filter,  yielding  an  effective  wavelength  of  4350  A. 
The  near  infrared  images  were  taken  using  hypersensitized  IV-N  emulsion  with  Wr  88 A 
filter,  yielding  an  effective  wavelength  of  8250  A.  Here  again,  the  reader  is  referred  to 
EE  for  further  details  concerning  the  observations  and  data  reduction. 

Primary  emphasis  has  been  placed  on  the  / passband  images  because  they  better 
describe  the  distribution  of  the  older  disk  stars,  the  gravitationally  dominant  component 
of  the  disk.  EE  determine  the  near-infrared  disk  scale  length  of  NGC  3992  to  be  3.38 
± 0.52  kpc,  assuming  a distance  of  11.3  Mpc  to  the  galaxy.  Converting  this  scale 
length  to  our  assumed  distance  of  14.2  Mpc  yields  a value  of  4.25  + 0.65  kpc,  or  an 
inverse  disk  scale  length  of  0.235 kpc  ^ EE  also  determined  the  slope  of  the 
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radial  dependence  of  the  spiral  arm-interarm  surface  brightness  difference,  from  which 
we  derive  the  inverse  scale  length  of  the  spiral  perturbation  of  0.4  kpc“  . 

Hunter  et  al.  (1988)  have  analyzed  the  bar  figure  of  NGC  3992  in  the  infrared 
image  taken  by  EE.  They  determined  the  orientation  and  apparent  extent  of  the  bar  from 
the  positions  of  the  extreme  ends  of  its  apparent  major  axis.  Further,  the  apparent  axial 
ratio  /5o  was  estimated  by  fitting  ellipses  to  the  outermost  isophotes  of  the  bar,  after 
subtraction  of  the  bulge  component.  They  found  that  the  major  axis  of  the  projected 
bar  figure  makes  an  angle  = 35°  with  respect  to  the  line  of  nodes  of  the  disk,  that  the 
semimajor  axis  of  the  bar  is  62",  and  that  the  apparent  axial  ratio  /?o  = 2.6.  These  three 
quantities  were  used  along  with  the  inclination  angle  i and  the  adopted  distance  to  derive, 
via  the  method  of  Stark  (1977,  also  see  Chapter  2),  the  actual  linear  dimensions  of  the 
triaxial  bar  figure  (again,  assuming  da  = 0.1).  The  resulting  long  bar  axis  a is  5.5  kpc, 
the  intermediate  axis  6 is  2.1  kpc,  and  the  short  axis  c is  0.55  kpc.  Again,  it  is  assumed 
that  the  long  and  intermediate  bar  axes  lie  in  the  disk  plane.  Finally,  Kennicutt  (1981) 
gives  the  spiral  arm  pitch  angle  of  NGC  3992  as  -11°  ± 2°.  Table  3-3  summarizes  the 
quantities  derived  from  the  surface  photometry  that  are  used  in  this  work. 

Best  Model 

The  first  step  in  determining  the  most  successful  model  of  NGC  3992  was  to  fit 
visually  the  observed  rotation  curve  with  the  theoretical  curve  comprised  of  contributions 
from  the  bar,  disk,  and  halo.  Next,  the  pattern  speed  ftp  was  fixed  by  the  assumption 
that  the  long  bar  axis  is  equal  to  the  corotation  radius.  Finally,  initial  values  of  the 
parameters  characterizing  the  bar/spiral  perturbation  were  set  by  the  surface  photometry. 
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Table  3-3:  Photometrically  derived  parameter  values  for  NGC  3992. 


Parameter 

Value 

Inverse  disk  scale  length  (kpc“')® 

n OQK  +0-043 

Inverse  spiral  scale  length  (kpc“')^ 

0.4 

Angle  between  projected  bar  major  axis  and  disk  line 
of  nodes  ('’)'’ 

35 

Apparent  bar  semimajor  axis  (arcsec)'’ 

62 

Apparent  bar  axial  ratio'’ 

2.6 

Lxjng  bar  axis  length  (kpc)'’ 

5.5 

Intermediate  bar  axis  length  (kpc)'’ 

2.1 

Short  bar  axis  length  (kpc)*’ 

0.55 

Spiral  arm  pitch  angle  (°)‘' 

-11  ± 2 

“Elmegreen  and  Elmegreen  (1985) 

'’Hunter,  Ball,  Huntley,  England,  and  Gottesman  (1988) 

‘’derived  using  the  method  of  Stark  (1977) 

‘'Kennicutt  (1981) 

The  periodic  orbits  of  the  model  were  then  determined,  and  the  surface  response  density 
was  calculated  and  measured  for  self-consistency  with  the  imposed  surface  density.  The 
free  model  parameters,  primarily  the  ones  characterizing  the  spiral  perturbation,  were 
adjusted  and  the  model  differentially  corrected  until  the  most  successful  model  was 
found.  Table  3-4  summarizes  the  parameters  of  the  most  successful  model  of  NGC 
3992. 

The  rotation  curve  of  this  model  is  shown  in  Figure  3-4.  Also  shown  are  the 
contributions  of  the  separate  model  components.  Figure  3-5  shows  the  characteristics 
of  the  different  orbit  families  which  are  included  in  the  model.  The  most  important 
periodic  families  in  the  bar  are  the  2/1  (Figure  3-6)  and  4/1  (Figure  3-7)  resonant 
families.  The  3/1  resonant  family  exists  over  a shorter  range  of  energy,  and  these  orbits 
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are  found  not  to  be  very  critical  to  the  overall  structure  of  the  bar.  Also  included  inside 
the  bar  radius  <3  is  a set  of  “circular”  orbits  extending  from  Vc  = 3.5  kpc  to  = 5.5  kpc. 


Table  3-A:  Parameter  values  for  the  best  self-consistent  model  of  NGC  3992. 


Ferrers  Bar 

Spiral 

Mb  = 1.5  X lO^'^M© 

A = 2000  km^  s"^  kpc“^ 

a = 5.5  kpc 

Ar  = 0 km^  s“^  kpc“^ 

b = 2.1  kpc 

Ss  = 0.4  kpc~^,  iQ  = -10° 

c = 0.55  kpc 

ri  = 1.5  kpc,  K2  = 10.6  kpc 

Vlp  = 43.6  km  s“^  kpc~^ 

= K2  = 1,  A = 0.1  kpc 

Exponential  Disk 

Plummer  Sphere  Halo 

CO  = 750M©  pc-2  Mh  = 2.75  x 10“M© 

£0  = 0.235  kpc-i  bH  = 12  kpc 


Figure  3-5:  Characteristics  of  the  orbit  families  included  in  the  model. 
Each  characteristic  plots  x,  where  a given  orbit  crosses  the  minor  bar 
axis  b,  as  a function  of  Jacobi  constant,  as  parameterized  by 
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Outside  the  bar  the  main  family  is  the  -2/1  resonant  family  (Figure  3-8).  These 
orbits  support  the  imposed  spiral  in  the  region  between  the  -4/1  resonance  (outside 
of  corotation)  and  the  outer  Lindblad  resonance  (OLR).  The  -4/1  and  -6/1  resonant 
families  are  also  included,  but  they  play  a much  less  significant  role  than  does  the 
-2/1  family  in  supporting  the  imposed  spiral  structure.  Two  more  sets  of  “circular” 
orbits  are  included  outside  of  corotation.  The  first,  a continuation  of  the  set  inside 
of  corotation,  extends  from  = 5.5  kpc  to  7.5  kpc.  The  second,  included  to  sample 
the  response  of  matter  near  and  slightly  beyond  the  OLR,  extends  from  rc  = 10  kpc 
to  12  kpc.  The  use  of  the  “circular”  orbits  near  corotation  does  not  imply  that  there 
are  no  stable  periodic  orbits  in  this  region.  However,  as  we  shall  see  in  Chapter  5, 
the  available  phase  space  in  this  interval  of  energy  is  dominated  by  stochasticity  and 
the  trapping  done  by  periodic  orbits  is  small.  Therefore,  given  the  magnitude  of  the 
imposed  velocity  dispersion,  it  makes  essentially  no  difference  to  the  resulting  density 
response  map  whether  these  orbits  have  the  initial  conditions  of  periodic  orbits  or  not. 
Furthermore,  despite  the  fact  that  the  majority  of  these  orbits  are  stochastic,  they  do 
provide  significant  enhancement  of  the  imposed  spiral  structure  and  are  instrumental  in 
achieving  self-consistency.  The  behavior  of  these  stochastic  orbits  will  be  considered 
in  more  detail  in  Chapter  5.  Table  3-5  summarizes  the  positions  of  the  main  resonances 
in  the  model.  The  epicyclic  frequency  «(r)  and  angular  rotation  rate  0(r),  necessary  for 
determining  the  resonance  positions,  are  obtained  by  averaging  and  - in 

azimuth,  where  Vsir,  9)  is  the  potential  of  the  bar,  and  adding  them  to  the  corresponding 
quantities  of  the  disk  and  halo.  k{v),  then,  is  simply  0(r)  is 

i 

, where  Vr(r)  is  the  total  axisymmetric  potential. 
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Figure  3-6:  The  2/1  family  of  periodic  orbits  in  the  model  of 
NGC  3992.  The  darker  circle  represents  corotation  at  5.5  kpc. 


Figure  3-7:  The  4/1  family  of  periodic  orbits  in  the  model  of  NGC  3992. 
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Figure  3-8:  The  -2/1  family  of  periodic  orbits  in  the  model  of  NGC  3992. 
The  darker  curves  represent  the  minima  of  the  bar  and  spiral  potentials. 


The  density  response  map  was  generated  according  to  the  procedure  outlined 
in  Chapter  2.  Table  3-6  gives  the  parameters  used  for  the  surface  density  response 
calculation. 


Table  3-5:  Resonance  locations  of  the  best  model  of  NGC  3992. 


Inner  Resonances 

Outer  Resonances 

Type 

Location 

Type 

Location 

2/1 

0.0  kpc 

-8/1 

6.9  kpc 

4/1 

3.1  kpc 

-6/1 

7.3  kpc 

6/1 

3.9  kpc 

-4/1 

8.2  kpc 

8/1 

4.3  kpc 

-2/1 

10.6  kpc 

corotation  = 5.5  kpc 
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Table  3-6:  Parameters  used  to  calculate  the  surface 
density  response  for  the  best  model  of  NGC  3992. 


Parameter 


Value 


Minimum  grid  radius,  RMIN  (kpc) 

Maximum  grid  radius,  RMAX  (kpc) 

Radial  cell  width,  DRS  (kpc) 

Number  of  radial  cells 
Azimuthal  cell  width  (°) 

Number  of  azimuthal  cells 

Number  of  radial  start  positions  per  grid  cell,  NOR 
Number  of  azimuthal  start  positions  per  orbit,  NOA 

Length  of  (quasi)periodic  orbit  integration,  PN 
Central  velocity  dispersion,  ao  (km  s“^) 

1 1 

Slope  of  the  velocity  dispersion  dependence,  Cr  (km  s“  kpc“  ) 
Time  of  "circular"  orbit  integration  (10^  yr) 


0.1 

12.1 

0.2 

60 

1.4 

128 

2 

6 

12 

100 

-7 

0.98 


Grayscale  representations  of  the  resulting  surface  density  response  are  shown  in 
Figures  3-9  (unprojected)  and  3-10  (projected  to  the  actual  orientation  of  NGC  3992,  for 
comparison  with  Figure  3-1).  These  grayscale  images  were  produced  by  first  rebinning 
the  polar  grid  of  density  values  into  a Cartesian  grid.  The  resulting  image  was  then 
smoothed  using  a “beam”  five  pixels  square  such  that  the  smoothed  pixel  derived  one- 
third  of  its  value  from  the  unsmoothed  pixel,  another  one-third  from  the  eight  directly 
adjacent  pixels,  and  the  final  third  from  the  sixteen  outer  pixels  (giving  a standard 
deviation  for  an  equivalent  Gaussian  beam  of  1.2  pixels).  There  are  a couple  of  features 
to  note  in  these  images.  Most  obviously  we  see  that  the  overall  barred  spiral  appearance 


of  NGC  3992  is  well-reproduced  by  the  model.  This  point  is  perhaps  best  appreciated  by 


directly  comparing  the  projected  figure  (Figure  3-10)  with  the  photograph  of  NGC  3992 
(Figure  3-1).  Also,  we  see  that  model  spiral  arm  makes  a slight  bend  at  approximately 
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the  position  of  the  -4/1  resonance  (8.2  kpc).  This  feature  seems  also  to  be  present 
in  NGC  3992  itself,  particularly  for  the  spiral  arm  in  the  northeast  quadrant.  These 
“elbows”  were  also  noted  by  Patsis  et  al.  (1991)  in  their  models  of  unbarred  spirals. 

After  the  surface  response  density  was  calculated,  the  self-consistency  of  the  model 
was  tested  by  comparing  the  response  to  the  imposed  surface  density.  Figures  3-11  and 
3-12  give  R and  A^,  the  two  measures  of  self-consistency.  These  measures,  the  ratio 
of  the  relative  amplitude  of  the  response  29  component  to  the  relative  amplitude  of 
the  imposed  29  component  and  the  phase  difference  between  these  two  components, 
respectively,  are  described  in  detail  in  Chapter  2. 

We  note  that  the  value  of  R is  almost  constant  at  one  (varying  between  1.4  and 
0.8)  between  2.1  and  8.5  kpc.  There  is  an  upward  divergence  in  R*  at  small  radii  (not 
shown  in  Figure  3-1 1)  due  to  the  fact  that  the  relative  azimuthal  variation  of  the  imposed 
density  drops  rather  rapidly  inside  the  length  of  the  bar  minor  axis,  greatly  decreasing  the 
relative  amplitude  of  the  imposed  29  component,  while  the  shape,  and  hence  the  density 
response,  of  the  dominant  2/1  family  remains  rather  elongated,  thereby  maintaining  the 
relative  amplitude  of  the  response  29  component.  These  combined  effects,  together 
with  the  possible  existence  of  a nuclear  bulge,  indicate  that,  in  general,  we  cannot 
meaningfully  speak  of  a measureable  bar  component  existing  all  the  way  to  the  center. 
Here,  therefore,  and  in  the  cases  of  NGC  1073  and  NGC  1398  also,  we  take  the  bar 
semiminor  axis  length  as  a practical  inner  limit  for  R . The  divergence  near  the  OLR 
is  due  to  several  factors.  First,  the  orbits  here  retain  their  somewhat  elongated  shape 
despite  the  fact  that  the  strength  of  the  imposed  spiral  has  considerably  diminished. 
Secondly,  there  exists  a congestion  of  matter  resulting  from  the  response  of  the  regular 
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Figure  3-9:  Grayscale  image  of  the  unprojected  surface 
density  response  of  the  best  model  of  NGC  3992. 
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Figure  3-10:  Grayscale  image  of  the  surface  density  response  of  the 
best  model  of  NGC  3992  projected  to  the  galaxy’s  actual  orientation. 
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orbits  trapped  (with  only  a small  dispersion  of  velocities)  around  the  -2/1  family. 
Thirdly,  and  perhaps  most  importantly,  the  effect  on  R of  an  azimuthal  density  variation 
of  a given  absolute  size  will  be  greatly  magnified  in  the  outer  regions  where  the 
axisymmetric  component  of  the  density  becomes  small. 

The  angle  deviation  A9  is  absolutely  less  than  0.12  radian  between  2.1  and  7.7 
kpc.  Beyond  7.7  kpc  we  see  that  the  response  spiral  systematically  lags  the  imposed  up 
to  the  -4/1  resonance.  Here  we  see  a short  sharp  increase  in  Ad  and  then  a much  steeper 
decline  out  to  the  OLR.  This  same  systematic  lag  has  been  observed  in  similar  models 
of  normal  spirals  with  strong  spiral  arms  (Contopoulos  and  Grosbpl  1988;  Patsis  et  al. 
1991),  and  seems  to  indicate  that  logarithmic  spirals  may  not  be  the  best  mathematical 
description  of  real  spirals. 


Figure  3-11:  The  response-to-imposed  29  component  amplitude 
ratio  R*  for  the  best  model  of  NGC  3992.  The  positions  of  the  major 
outer  resonances  (corotation,  -A!\ , and  outer  Lindblad)  are  noted. 
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kpc 

Figure  3-12:  The  phase  difference  A9  (in  radians)  between  the  response 
and  imposed  20  components  of  the  best  model  of  NGC  3992.  The 
positions  of  the  major  outer  resonances  are  noted  as  in  Figure  3-11. 


NGC  1073 


Observations 


Neutral  Hydrogen 


In  a continuation  of  the  effort  initiated  by  GBHH  to  study  the  properties  of 
neutral  hydrogen  in  a carefully  selected  sample  of  barred  spiral  galaxies,  NGC  1073 
was  observed  by  England,  Gottesman,  and  Hunter  (1990,  hereafter  EGH)  at  the  VLA 
between  June,  1983  and  June,  1984.  The  effective  resolution  obtained  was  20. "3  x 19. "7, 
with  the  major  axis  of  the  Gaussian  beam  at  a position  angle  of  62.°  8.  Eor  comparison 
purposes,  the  diameter  of  the  HI  disk  is  396"  and  the  position  angle  of  the  line  of  nodes 
is  -15. °4  ± 0.°2.  The  FWHP  velocity  resolution  obtained  was  12.63  km  s~^.  EGH 
gives  full  details  of  the  VLA  data  reduction. 
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The  angle-averaged  rotation  curve,  along  with  the  following  orientation  parameters 
of  the  HI  disk,  was  derived  via  a least-squares  fit  to  the  velocity  field:  major  axis  at 
a position  angle  of-15.°4  + 0.°2,  inclination  of  18.°5  ± 2.°5,  and  a systemic  velocity 
of  1208.9  + 0.2  km  s"^.  The  rotation  curve  data  is  given  in  Table  3-7.  Here  again, 
the  quoted  uncertainties  are  simply  the  formal  errors  of  the  least-squares  solution.  In 
converting  the  measured  angular  radii  to  a linear  scale  we  have  assumed  a distance  to 
NGC  1073  of  13.6  Mpc  (de  Vaucouleurs  1979).  The  rotation  curve  data  in  Table  3-7 
are  plotted  as  filled  circles  with  error  bars  in  Figure  3-13. 


Table  3-7:  HI  rotation  velocity  data  for  NGC  1073. 


R (arcmin) 
(1) 

R (kpc) 
(1) 

V (km  s-^) 
(2) 

R (arcmin) 
(1) 

R (kpc) 
(1) 

V (km  s ^) 
(2) 

0.10 

0.39 

38.8  ± 7.6 

2.70 

10.6 

100.8  + 5.6 

0.30 

1.18 

83.2  ± 4.4 

2.90 

11.4 

95.2  ± 6.6 

0.50 

1.97 

96.0  + 3.2 

3.10 

12.2 

92.0  + 7.7 

0.70 

2.76 

93.6  ± 2.4 

3.30 

13.0 

100.8  + 8.4 

0.90 

3.55 

98.8  + 3.2 

3.50 

13.8 

91.2  + 12.0 

1.10 

4.34 

99.6  ± 3.2 

3.70 

14.6 

69.6  ± 26.4 

1.30 

5.12 

99.6  ± 2.2 

3.90 

15.4 

74.4  + 28.3 

1.50 

5.91 

100.0  ± 2.8 

4.10 

16.2 

61.6  + 31.6 

1.70 

6.70 

105.6  ± 3.2 

4.30 

17.0 

24.0  + 33.0 

1.90 

7.49 

106.0  + 3.2 

4.50 

17.7 

72.0  + 44.4 

2.10 

8.28 

103.2  ± 2.8 

4.70 

18.5 

120.0  ± 52.4 

2.30 

9.07 

104.0  + 4.4 

4.90 

19.3 

168.0  ± 53.6 

2.50 

9.85 

105.6  ± 6.0 

% 


(1)  Radius  measured  from  the  center  of  the  galaxy 

(2)  Rotation  velocity 
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EGH  find  that  the  distribution  of  the  neutral  hydrogen  in  NGC  1073  is  that  of 
an  almost  circular  disk,  with  a nearly  complete  ring  surrounding  a central  depression; 
however,  this  central  depression  is  not  as  deep  relative  to  the  outer  disk  as  in  the  case 
of  NGC  3992.  They  also  note  the  presence  of  a gas  bar  in  NGC  1073  aligned  with  the 
optical  bar.  This  feature  shows  up  at  around  30%  of  the  peak  HI  density,  compared 
to  the  general  disk  gas  density  which  is  approximately  20%  of  peak.  EGH  find  little 
evidence  of  spiral  structure  in  the  gas,  however,  in  sharp  contrast  to  the  optical  image 
in  which  there  are  two  prominent  spiral  arms  starting  about  30°  in  azimuth  from  the 
ends  of  the  bar.  While  the  regions  of  high  HI  density  (>50%  of  the  peak  density)  show 
a correlation  with  the  bright  optical  regions  of  the  spiral  arms,  a broad  gaseous  ring 
emerges  when  lower  HI  densities  (-40%  of  peak)  are  considered.  This  ring  extends 
from  about  r = 0.'7  to  r = 2.'0,  compared  to  the  optical  spiral  arms  which  extend  from 
about  r = 0.'75  to  r = l.'l.  Finally,  EGH  note  a steep  gradient  in  the  gas  density  in  the 
northwest  quadrant  of  the  galaxy.  While  they  postulate  that  this  feature  might  be  due 
to  an  interaction,  a search  for  21  cm  emission  in  a l.°5  x l.°5  area  of  sky  surrounding 
NGC  1073  yielded  no  evidence  of  such  an  interacting  object. 

The  gas  kinematics  of  NGC  1073  is  dominated  by  circular  motion.  Still,  irregular- 
ities in  the  velocity  field  due  to  gas  streaming  are  observed  where  the  optical  spiral  arms 
cross  the  isovelocity  contours  (EGH  1990).  Like  NGC  3992,  NGC  1073  may  possess  a 
truncated  disk.  EGH  cite  as  evidence  for  this  the  dropoff  in  the  rotation  curve  starting 
at  approximately  r = 10  kpc,  the  corresponding  drop  in  the  angle-averaged,  deprojected 
HI  surface  density,  and  the  success  achieved  by  truncated  mass  models  in  reproducing 
the  rotation  curve  of  NGC  1073  (cf.  Casertano  1983;  Hunter  et  al.  1984). 
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Figure  3-13;  Comparison  of  observed  and  theoretical  rotation  curves  for  NGC  1073. 


Table  3-8:  Selected  global  and  disk  parameter  values  adopted  for  NGC  1073.  Here  again,  the 
errors  listed  are  simply  the  formal  errors  of  a least-squares  analysis  of  the  velocity  field. 


Parameter 

Value 

Systemic  velocity  (heliocentric,  km  s“^)^ 

1208.9  + 0.2 

Distance  (Mpc)’’ 

13.6 

Inclination  angle  (°Y 

18.5  + 2.5 

Position  angle,  line  of  nodes  (°)^ 

-15.4  + 0.2 

^England,  Gottesman,  and  Hunter  (1990) 
’’de  Vaucouleurs  (1979) 

Surface  Photometry 

As  for  NGC  3992,  Elmegreen  and  Elmegreen  (1985)  have  provided  B and  / 
passband  surface  photometry  of  NGC  1073.  They  determined  the  near-infrared  disk 
scale  length  of  NGC  1073  to  be  3.18  ± 0.36  kpc  assuming  a distance  of  13.2  Mpc.  For 
our  assumed  distance  of  13.6  Mpc,  the  scale  length  becomes  3.28  ± 0.37  kpc,  yielding 
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an  inverse  disk  scale  length  of  0.305 ^Q  ggj  kpc“^  The  initial  value  of  the  inverse  scale 
length  of  the  spiral,  derived  from  the  slope  of  the  spiral  arm-interarm  surface  brightness 
difference  given  by  EE,  is  1.4  kpc“^ 

EGH  have  analyzed  the  bar  figure  of  NGC  1073  in  the  infrared  image  taken  by 
EE.  Following  the  same  procedure  described  for  the  case  of  NGC  3992,  they  find  that 
the  angle  ^ between  the  major  axis  of  the  projected  bar  and  the  line  of  nodes  is  75°. 
Also,  they  find  that  the  bar  has  a semimajor  axis  = 43"  and  an  apparent  axial  ratio  j3o 
= 7.19.  The  actual  linear  dimensions  of  the  bar  derived  from  the  Stark  method  are  a 
= 2.95  kpc,  b = 0.39  kpc,  and  c = 0.295  kpc.  The  initial  value  of  the  spiral  arm  pitch 
angle,  determined  by  a visual  fit  to  the  optical  photograph  of  NGC  1073  reproduced 
in  EGH  from  Arp  and  Sulentic  (1979),  is  -10°.  Table  3-9  summarizes  the  quantities 
derived  from  the  surface  photometry. 


Table  3-9:  Photometrically  derived  parameter  values  for  NGC  1073. 


Parameter 

Value 

Inverse  disk  scale  length  (kpc~')“ 

n QnK+O-039 
U.OUO_Q  Q3^ 

Inverse  spiral  scale  length  (kpc“')^ 

1.4 

Angle  between  projected  bar  major  axis  and  disk  line 
of  nodes  (°)'’ 

75 

Apparent  bar  semimajor  axis  (arcsec)'’ 

43 

Apparent  bar  axial  ratio’’ 

7.19 

Long  bar  axis  length  (kpc)° 

2.95 

Intermediate  bar  axis  length  (kpc)° 

0.39 

Short  bar  axis  length  (kpc)° 

0.295 

Spiral  arm  pitch  angle  (°)*' 

-10 

^Elmegreen  and  Elmegreen  (1985) 

'’England,  Gottesman,  and  Hunter  (1990) 
^derived  using  method  of  Stark  (1977) 

‘'visual  fit 
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Best  Model 


The  best  model  for  NGC  1073  was  determined  following  the  general  procedure 
outlined  for  the  case  of  NGC  3992.  Table  3-10  summarizes  the  parameters  of  this  model. 


Table  3—10:  Parameter  values  for  the  best  self-consistent  model  of  NGC  1073. 


Ferrers  Bar 

Ms  = 1.8  X lO^M© 
a = 2.95  kpc 

b = 0.39  kpc 

c = 0.295  kpc 

ftp  = 32.2  km  s“^  kpc“^ 


Spiral 

A = 9000  km^  s”^  kpc"^ 

Ar  = 0 km^  s“^  kpc"^ 

= 1.2  kpc“^  io  = -10° 
ri  = 2.95  kpc,  V2  = 5.63  kpc 
Ki  = K2  = 0.5,  A = 0.1  kpc 


Exponential  Disk  Plummer  Sphere  Halo 

CO  = 250M©  pc"2  Mh  = 1.0  X lO^^M© 

£0  = 0.305  kpc-i  bn  = 9 kpc 

The  rotation  curve  of  this  model  is  given  along  with  the  observations  in  Figure 
3-13.  The  characteristics  of  the  orbit  families  which  comprise  this  model  are  shown  in 
Figure  3-14.  The  most  dynamically  important  periodic  orbits,  that  is,  those  that  are  most 
effective  in  trapping  quasi-periodic  orbits  around  them,  are  the  two  stable  branches  of 
the  2/1  resonant  orbits  in  the  bar  and  the  -2/1  resonant  orbits  in  the  outer  disk  (Figure 
3-15).  Two  sets  of  “circular”  orbits  are  also  included,  a sizeable  one  extending  from 
Kc  = 1.25  kpc  to  Tc  = 4.75  kpc  and  a smaller  second  one  extending  from  Vc  = 5.25  kpc 
to  Vc  = 5.6  kpc.  As  in  the  case  of  NGC  3992,  the  inclusion  of  these  “circular”  orbits 
does  not  imply  that  no  periodic  orbits  exist  in  this  region.  We  shall  explore  further  the 
phase  space  structure  of  this  region  in  Chapter  5.  Table  3-11  gives  the  positions  of  the 


main  resonances  in  this  case. 


78 


Figure  3-14:  Characteristics  of  the  orbit  families  included  in  the  best  model  of  NGC  1073. 


kpc 

Figure  3-15:  Representative  orbits  of  the  three  main 
periodic  families  in  the  best  model  of  NGC  1073. 
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Table  3-11:  Resonance  locations  of  the  best  model  of  NGC  1073. 


Inner  Resonances 

Outer  Resonances 

Type 

Location 

Type 

Location 

2/1 

0.0  kpc 

-8/1 

3.6  kpc 

4/1 

1.8  kpc 

-6/1 

3.8  kpc 

6/1 

2.2  kpc 

-4/1 

4.3  kpc 

8/1 

2.4  kpc 

-2/1 

5.6  kpc 

corotation  = 2.95  kpc 


Table  3-12  gives  the  parameters  of  the  surface  density  response  calculation. 


Table  3-12:  Parameters  used  to  calculate  the  surface 
density  response  for  the  best  model  of  NGC  1073. 

Parameter 

Minimum  grid  radius,  RMIN  (kpc) 

Maximum  grid  radius,  RMAX  (kpc) 

Radial  cell  width,  DRS  (kpc) 

Number  of  radial  cells 
Azimuthal  cell  width  (°) 

Number  of  azimuthal  cells 

Number  of  radial  start  positions  per  grid  cell,  NOR 
Number  of  azimuthal  start  positions  per  orbit,  NOA 

Length  of  (quasi)periodic  orbit  integration,  PN  f 2*x  N oA  ) 

Central  velocity  dispersion,  ao  (km  s“^) 

Slope  of  the  velocity  dispersion  dependence,  <t^  (km  s“^  kpc“^) 
Time  of  "circular"  orbit  integration  (10^  yr) 


Value 

0.02 

6.02 

0.1 

60 

1.4 

128 

2 

6 

12 

40 

-5 

0.98 


An  unprojected  grayscale  representation  of  the  surface  density  response  is  given 
in  Figure  3-16.  Since  the  inclination  of  the  galaxy  is  only  18.°5  (i.e.  it  is  nearly  face 
on),  a projected  image  is  not  given  here.  As  in  the  case  of  NGC  3992,  we  see  that  the 
overall  barred  spiral  structure  is  well-reproduced.  Also,  in  contrast  to  NGC  3992  and  its 
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Figure  3-16:  Grayscale  image  of  the  unprojected  surface 
density  response  of  the  best  model  of  NGC  1073. 
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best  model,  the  spiral  arms  of  this  model  do  not  exhibit  any  noticeable  “elbows.” 
Furthermore,  it  is  observed  that  the  model  arms  tend  to  wrap  back  onto  themselves  after 
winding  approximately  270°  from  their  origin  at  the  ends  of  the  bar.  This  “wrapping” 
occurs  at  about  the  radius  of  the  OLR.  One  note  of  discrepancy  between  the  model  and 
observations  exists,  though.  While  the  actual  spiral  arms  of  NGC  1073  are  observed  to 
begin  some  30°  in  azimuth  from  the  ends  of  the  bar  (EGH  1990),  the  arms  of  the  best 
model  come  directly  off  the  ends  of  the  bar. 


Figure  3-17:  The  response -to-imposed  20  component  amplitude 
ratio  R*  for  the  best  model  of  NGC  1073.  The  positions  of  the  major 
outer  resonances  (corotation,  -4/1 , and  outer  Lindblad)  are  noted. 


Figures  3-17  and  3-18  give  R*  and  AO  for  the  best  model  of  NGC  1073.  Again 
we  see  the  same  general  features  here  as  in  the  case  of  NGC  3992.  The  value  of  R*  is 
almost  constant  (varying  between  0.8  and  1.3)  from  0.6  kpc  to  4.1  kpc.  There  is  also, 
for  the  same  reasons  as  in  the  model  of  NGC  3992,  an  upward  divergence  in  R*  at  large 
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radii.  The  angle  deviation  /S.6  remains  absolutely  less  than  0.12  radian  from  0.39  kpc 
out  to  3.9  kpc.  Again  we  see  a gradual  systematic  lag  in  the  response  26  component 
relative  to  the  imposed  from  the  end  of  the  bar  at  corotation  out  to  approximately  the 
-4/1  resonance.  Beyond  this  point  A0  drops  more  rapidly  out  to  the  OLR,  albeit  not 
as  sharply  as  in  the  case  of  NGC  3992. 


kpc 

Figure  3-18:  The  phase  difference  (in  radians)  between  the 
response  and  imposed  20  components  of  the  best  model  of  NGC 
1073.  Again,  the  positions  of  the  major  outer  resonances  are  noted. 


NGC  1398 


Observations 


Neutral  Hydrogen 


Moore  and  Gottesman  (1993)  made  high  resolution  HI  observations  of  NGC  1398 
at  the  VLA  between  February,  1991  and  February,  1992.  The  effective  linear  resolution 
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obtained  was  22."3  x 20. "2,  with  the  major  axis  of  the  Gaussian  beam  at  a position 
angle  of  -35.°9.  The  diameter  of  the  HI  disk  was  found  to  be  approximately  720"  and 
the  position  angle  of  the  line  of  nodes  is  -84.°  5 + 0.°5.  The  FWHP  velocity  resolution 
was  25.4  km  s“^. 

From  the  neutral  hydrogen  velocity  maps  an  angle-averaged  rotation  curve  for 
NGC  1398  was  derived  in  the  same  way  as  for  NGC  3992  and  NGC  1073.  In  the 
fitting  procedure  the  HI  disk  was  determined  to  have  the  following  properties:  major 
axis  at  a position  angle  of  -84.°5  ± 0.°5,  inclination  of  46.°3  ± 0.°7,  and  a systemic 
velocity  of  1397.0  ±1.0  km  s“^.  As  in  the  cases  of  NGC  3992  and  NGC  1073,  the  HI 
surface  density  distribution  of  NGC  1398  shows  evidence  of  a depletion  in  the  central 
region.  Moore  (private  communication)  notes  that  this  HI  “hole”  is  rather  small  (at  a 
lower  spatial  resolution  — 32"  x 60"  beam)  extending  about  30"  to  40"  from  center,  a 
distance  only  slightly  less  than  the  deprojected  semimajor  axis  of  the  bar.  Also,  in  the 
high  resolution  HI  density  map,  a rather  clumpy  gas  ring  extending  over  the  region  of 
the  outer  optical  spiral  arms  is  observed.  Analysis  of  the  high  resolution  data  yielded  a 
rotation  curve  which  contained  no  information  about  the  central  rising  part  of  the  curve. 
Therefore,  another  rotation  curve  was  then  derived  using  the  lower  spatial  resolution 
in  order  to  increase  the  signal-to-noise  ratio  in  the  inner  region.  The  rotation  curve 
finally  adopted  for  NGC  1398  is,  then,  a composite  of  the  two.  In  the  inner  region  the 
rotation  data  come  from  analysis  of  the  low  resolution  map  only.  For  the  remainder 
of  the  rotation  curve  an  average  of  the  two  data  sets  is  used.  Near  the  edge  of  the  HI 
disk  the  two  curves  diverge  for  reasons  unknown,  but  this  occurs  beyond  the  end  of 
our  adopted  rotation  curve  and  does  not  affect  our  results.  The  rotation  curve  data  are 
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given  in  Table  3-13.  The  distance  to  NGC  1398  is  taken  to  be  18.7  Mpc  (Buta  and  de 
Vaucouleurs  1983).  The  data  are  plotted  as  filled  circles  with  error  bars  in  Figure  3-19. 


Table  3-13:  HI  rotation  velocity  data  for  NGC  1398. 


R (arcmin) 
(1) 

R (kpc) 
(1) 

V (km  s-^) 
(2) 

R (arcmin) 
(1) 

R (kpc) 
(1) 

V (km  s ^) 
(2) 

0.26 

1.41 

128.2  + 65.7 

2.97 

16.2 

301.5  + 5.0 

0.54 

2.94 

250.6  + 19.8 

3.32 

18.1 

295.4  ± 5.4 

0.89 

4.84 

293.3  + 5.8 

3.67 

20.0 

291.2  + 6.4 

1.24 

6.75 

299.1  ± 4.9 

4.03 

21.9 

289.2  ± 4.1 

1.59 

8.65 

300.7  ± 4.3 

4.37 

23.8 

289.6  ± 1.3 

1.93 

10.5 

307.7  ± 4.6 

4.72 

25.7 

285.9  ± 2.5 

2.28 

12.4 

307.7  + 2.9 

5.06 

27.5 

279.8  + 6.4 

2.64 

14.4 

305.2  + 3.8 

(1)  Radius  measured  from  the  center  of  the  galaxy 

(2)  Rotation  velocity 


Figure  3-19:  Comparison  of  observed  and  theoretical  rotation  curves  for  NGC 
1398.  The  contributions  of  the  separate  model  components  are  also  shown. 
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Table  3-14;  Selected  global  and  disk  parameter  values  adopted  for  NGC  1398. 


Parameter 

Value 

Systemic  velocity  (heliocentric,  km  s“^)^ 

1397  + 1 

Distance  (Mpc)*’ 

18.7 

Inclination  angle  (°)'^’‘^ 

43 

Position  angle,  line  of  nodes  (°)*^^ 

99.8 

^Moore  and  Gottesman  (1993) 

'’Buta  and  de  Vaucouleurs  (1983) 

^Ohta,  Hamabe,  and  Wakamatsu  (1990);  this  source  quotes  no  errors  for  these 
quantities. 

^*These  values  are  used  since  the  corresponding  values  from  Moore  and  Gottesman 
(1993)  were  not  available  at  the  time  of  the  model  runs. 

Surface  Photometry 


Ohta  et  al.  (1990)  have  provided  B passband  photometry  of  six  early-type  barred 
spiral  galaxies,  including  NGC  1398.  They  have  taken  two  photographs  of  NGC  1398 
(one  6 minute  exposure  and  one  30  minute  exposure)  with  the  2.5  m du  Pont  telescope 
at  the  Las  Campanas  Observatory  using  hypersensitized  103a-O  emulsions.  The  reader 
is  referred  to  Ohta  et  al.  (1990)  for  complete  details  on  the  observations  and  data 
reduction. 

They  decomposed  the  azimuthally  averaged  luminosity  profile  of  NGC  1398  into 
an  exponential  disk  and  an  r^^'^-law  bulge  by  an  iterative  scheme.  They  determine  the 
scale  length  of  the  exponential  disk  to  be  59."6.  Converting  this  scale  length  to  linear 
dimensions  by  using  our  assumed  distance  yields  an  inverse  scale  length  of  0.185  kpc“^. 
This  is  the  value  we  have  adopted  here.  We  were  unable  to  find  published  data  regarding 
the  scale  length  of  the  spirals,  thus  we  left  the  inverse  spiral  scale  length  parameter  to 


be  fit  in  the  modeling  procedure. 
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They  also  analyzed  the  bar  figure  of  NGC  1398  in  their  study.  They  give 
inclination-corrected  values  for  the  half-length  and  full  width  of  the  bar  of  53"  and  18", 
respectively.  However,  NGC  1398  possesses  a significant  bulge  component  and  these 
bar  dimensions  are  derived  only  from  portions  of  the  bar  which  extend  beyond  the  bulge 
radius  (i.e.  where  the  light  contributions  from  the  two  components  are  approximately 
equal),  which  they  determine  to  be  22"  for  NGC  1398.  This  method  overestimates  the 
axial  ratio  of  the  bar  in  cases,  such  as  this  thesis,  where  only  one  component  is  used 
to  represent  the  bar.  Instead,  the  following  method  was  used  to  derive  the  observed 
axial  ratio  !3q.  Ohta  et  al.  (1990)  also  give  isophotes  of  the  extracted  bars  (i.e.  with 
the  axisymmetric  components  subtracted,  see  Ohta  et  al.  [1990]  for  details  of  the 
extraction  process).  The  published  isophotes  for  the  bar  of  NGC  1398  derived  in  this 
way  were  measured  for  axial  ratio,  and  an  average  value  of  = 2.15  was  obtained. 
The  given  deprojected  half-length  of  the  bar  was  then  reprojected,  yielding  a projected 
bar  semimajor  axis  of  39".  Also,  from  the  isophotes  of  the  extracted  bar  we  measured 
the  angle  between  the  projected  bar  and  the  line  of  nodes  to  be  i/’  = 91°.  These  three 
quantities,  along  with  the  adopted  distance  to  NGC  1398  and  the  assumption  that  da 
= 0.1,  were  then  used  in  the  Stark  analysis  to  derive  the  actual  linear  dimensions  of 
the  triaxial  model  bar.  Finally,  Kennicutt  (1981)  gives  the  spiral  arm  pitch  angle  of 
NGC  1398  as  -6°  ± 2°.  Table  3-15  summarizes  the  parameters  derived  for  NGC  1398 
from  the  surface  photometry. 

Best  Model 

The  best  model  of  NGC  1398  was  determined  by  the  same  method  as  the  best 
models  for  NGC  3992  and  NGC  1073.  Table  3-16  summarizes  its  parameters. 
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Table  3-15:  Photometrically  derived  parameter  values  for  NGC  1398. 


Parameter 

Value 

Inverse  disk  scale  length  (kpc“^)^ 

0.185 

Angle  between  projected  bar  major  axis  and  disk  line 
of  nodes  (°)^ 

91  ± 1 

Apparent  bar  semimajor  axis  (arcsec)® 

39 

Apparent  bar  axial  ratio^ 

2.15 

Long  bar  axis  length  (kpc)'’ 

4.8 

Intermediate  bar  axis  length  (kpc)'’ 

1.6 

Short  bar  axis  length  (kpc)*’ 

0.48 

Spiral  arm  pitch  angle  {°y 

— 6 ± 2 

®Ohta,  Hamabe,  and  Wakamatsu  (1990) 

^derived  using  method  of  Stark  (1977) 
‘^Kennicutt  (1981) 

Table  3-16:  Parameter  values  for  the  best  self-consistent  model  of  NGC  1398. 

Ferrers  Bar 

Spiral 

Mb  = 1.2  X 10^°M© 
a = 4.8  kpc 
b = 1.6  kpc 
c = 0.48  kpc 
Qp  = 54.2  km  s“^  kpc"^ 

A = 3000  km^  s“^  kpc“^ 

Ar  = 0 km^  s“^  kpc“^ 

€s  = 0.4  kpc"\  i‘o  = -6° 
r\  = 1.0  kpc,  f2  » OLR 
K\  = K2  = 0.5,  A = 0.1  kpc 

Exponential  Disk 

Plummer  Sphere  Halo 

Co  --  1475Mq  pc“^ 
£o  = 0.185  kpc“^ 

Mh  = 6.0  X 10“  M© 
bn  = 35  kpc 
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Figure  3-20:  Characteristics  of  the  orbit  families  included  in  the  best  model  of  NGC  1398. 

The  rotation  curve  of  this  model  is  given  along  with  the  observations  in  Figure 
3-19.  The  characteristics  of  the  orbit  families  which  comprise  this  model  are  shown 
in  Figure  3-20.  Once  again  the  most  important  periodic  family  in  the  bar  is  the  2/1 
resonant  family.  There  exists  in  the  model  a branch  of  4/1  resonant  orbits,  but  they 
are  of  less  dynamical  importance  than  the  2/1  family.  This  is  because  they  exist  over  a 
much  smaller  range  of  energy  and  trap  much  less  matter  around  them  than  do  the  2/1 
orbits.  Also,  while  there  is  a very  small  family  of  -2/1  orbits  at  about  Vc  = 12.0  kpc, 
the  dominant  orbital  behavior  in  the  outer  disk  of  this  model  is  decidedly  stochastic. 
Accordingly,  we  include  a set  “circular”  orbits  extending  all  the  way  from  = 3.3  kpc 
to  Vc  = 11.8  kpc.  We  shall  examine  this  stochasticity  further  in  Chapter  5.  Figure  3-21 
shows  representative  orbits  from  the  three  stable  periodic  families  included  in  the  model. 
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Figure  3-21:  Representative  orbits  of  the  three  main 
periodic  families  in  the  best  model  of  NGC  1398. 


Table  3-17  gives  the  positions  of  the  major  resonances  of  the  model. 


Table  3-17:  Resonance  locations  of  the  best  model  of  NGC  1398. 


Inner  Resonances 

Outer  Resonances 

Type 

1 Location 

Type 

Location 

2/1 

0.0  kpc 

-8/1 

6.1  kpc 

4/1 

2.5  kpc 

-6/1 

6.6  kpc 

6/1 

3.2  kpc 

-4/1 

7.5  kpc 

8/1 

3.6  kpc 

-2/1 

9.8  kpc 

corotation  = 4.8  kpc 


Table  3-18  summarizes  the  parameters  used  in  the  generation  of  the  surface  density 


response  of  the  model. 
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Table  3-18:  Parameters  used  to  calculate  the  surface 
density  response  for  the  best  model  of  NGC  1398. 

Parameter 

Minimum  grid  radius,  RMIN  (kpc) 

Maximum  grid  radius,  RMAX  (kpc) 

Radial  cell  width,  DRS  (kpc) 

Number  of  radial  cells 
Azimuthal  cell  width  (°) 

Number  of  azimuthal  cells 

Number  of  radial  start  positions  per  grid  cell,  NOR 
Number  of  azimuthal  start  positions  per  orbit,  NO  A 

Length  of  (quasi)periodic  orbit  integration,  PN  j 

Central  velocity  dispersion,  ao  (km  s“^) 

Slope  of  the  velocity  dispersion  dependence,  ar  (km  s“^  kpc“^) 
Time  of  "circular"  orbit  integration  (10^  yr) 


Value 

0.01 

12.01 

0.2 

60 

1.4 

128 

2 

6 

12 

80 

-5 

0.98 


Figure  3-22  shows  a grayscale  representation  of  the  surface  density  response  for 
the  best  model  of  NGC  1398.  Since  the  bar  is  almost  perpendicular  to  the  line  of 
nodes  and  the  inclination  angle  is  not  very  large  (thus  the  main  effect  of  the  galaxy’s 
inclination  is  merely  to  foreshorten  the  bar),  we  do  not  present  a projected  grayscale 
image.  As  Figure  3-22  shows,  the  model  exhibits  a well-defined  structure,  with  the 
outer  spirals  extending  smoothly  outwards  with  no  “elbows”  and  no  strong  tendency  to 
wrap  back  onto  themselves.  This  feature  is  consistent  with  the  observed  optical  arms 
(Figure  3-3)  which  are  long  and  continuous  and  “can  be  traced  for  about  one  and  a 
quarter  revolutions  from  their  origin  ...”  (Sandage  1961).  One  feature  of  NGC  1398 
which  is  not  well-reproduced  by  our  model  is  the  inner  ring,  composed  of  nearly  circular 
spiral  segments,  which  encircles  the  bar.  It  is  not  clear,  however,  that  this  is  actually  a 
significant  feature  in  the  distribution  of  the  underlying  disk  stars. 


91 


Figure  3-22:  Grayscale  image  of  the  unprojected  surface 
density  response  of  the  best  model  of  NGC  1398. 


10 
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Figures  3-23  and  3-24  give  R*  and  AO  for  the  most  successful  model  of  NGC 
1398.  In  this  case  we  do  not  have  such  a pronounced  divergence  of  R at  large  radii 
as  we  had  in  the  cases  of  NGC  3992  and  NGC  1073.  Except  for  a small  dip  to  0.4 
around  9.4  kpc,  R remains  fairly  constant  between  0.8  and  1.3  from  1.6  kpc  to  11.2 
kpc.  The  angle  deviation  AO  remains  absolutely  less  than  0.4  radian  all  the  way  from 
1.6  kpc  out  to  10.6  kpc.  While  this  value  is  somewhat  larger  than  in  the  cases  of  NGC 
3992  and  NGC  1073,  the  agreement  is  still  rather  good  because  the  angle  deviation  is 
magnified  by  the  small  absolute  value  of  the  pitch  angle  (i.e.  given  the  same  deviation 
perpendicular  to  the  spiral,  the  angular  difference  will  be  greater  at  a specified  radius 
the  lower  is  the  pitch  angle).  In  this  case  again,  we  see  the  tendency  for  the  response 
spiral  progessively  to  lag  the  imposed. 


Figure  3-23:  The  response-to-imposed  20  component  amplitude 
ratio  R*  for  the  best  model  of  NGC  1398.  The  positions  of  the  major 
outer  resonances  (corotation,  -A/\,  and  outer  Lindblad)  are  noted. 
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Figure  3-24:  The  phase  difference  A9  (in  radians)  between  the 
response  and  imposed  26  components  of  the  best  model  of  NGC 
1398.  Again,  the  positions  of  the  major  outer  resonances  are  noted. 


CHAPTER  4 

VARIATION  OF  PARAMETERS 


In  this  chapter  we  examine  how  systematically  varying  the  parameters  affects  the 
self-consistency  of  the  “best”  models  described  in  Chapter  3.  The  purpose  of  such  an 
examination  is  twofold.  Firstly,  we  would  like  to  identify  those  parameters  which  exert 
the  most  influence  on  self-consistency.  In  this  way  we  might  hope  to  identify  quantities 
which  act  as  eigenvalues  of  the  self-consistency  problem.  Secondly,  we  would  like  to 
detect,  if  possible,  any  systematic  trends  in  this  set  of  controlling  parameters  as  we 
move  from  early  to  late  Hubble  types.  Since  varying  the  same  parameter  in  different 
models  often  produces  essentially  the  same  effect,  we  present  complete  details  for  the 
case  of  NGC  3992,  and  include  the  results  of  NGC  1073  and  NGC  1398  only  where 
they  differ  significantly  from  those  of  NGC  3992. 

Variation  of  A 

Figure  4—1  shows  how  varying  the  amplitude  A of  the  bar/spiral  perturbation 
potential  affects  self-consistency.  ^ Decreasing  the  amplitude  leads  to  worse  self- 
consistency  as  regards  both  R*  (at  all  radii)  and  AO  (outside  corotation).  On  the  other 
hand,  while  increasing  A produces  somewhat  smoother  runs  of  R and  AO  with  radius, 

R is  seen  to  decrease  systematically  below  unity  outside  of  corotation. 

1.  In  each  of  the  figures  in  this  chapter  we  have  included,  for  the  sake  of  completeness,  the  mns 
of  R (or  A9)  versus  radius  for  the  entire  radial  range  spanned  by  the  computational  grid.  For 
reasons  detailed  in  Chapter  3,  though,  we  do  not  expect  good  agreement  over  this  entire  radial 
range.  Specifically,  the  upward  divergence  noted  at  very  small  radii  can  be  safely  ignored. 
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Figure  4-1:  The  ratio  R*  and  phase  difference  A0  in  models  of  NGC  3992  where  A = 2000 
km^s“^kpc~^  (“best”  model,  denoted  by  filled  circles),  A = 1000  km^s“^kpc“^  (open  circles), 
and  A = 4000  km^s“^kpc“^  (plusses). 
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Variation  of  i4^ 


The  only  galaxy  whose  “best”  model  shows  any  real  sensitivity  to  variations  in 
the  residual  amplitude  Ar  is  NGC  1073.  The  reason  is  that  in  the  models  of  NGC  3992 
and  NGC  1398  the  inner  and  outer  cutoff  radii  (ri  and  ^2)  for  the  primary  amplitude 
A are,  respectively,  close  to  the  center  and  at  or  beyond  the  OLR.  That  is,  in  these 
cases  the  imposed  outer  spiral  perturbation  is  continued  well  inside  of  corotation  as  a 
bar  perturbation  in  addition  to  the  Ferrers  bar  component.  Therefore,  as  the  primary 
amplitude  extends  over  most  of  the  computational  grid  in  these  models,  the  residual 
amplitude  does  not  really  come  into  play.  The  “best”  model  of  NGC  1073,  however, 
has  the  inner  cutoff  r\  located  at  the  corotation  resonance  (i.e.  the  end  of  the  bar),  and 
variations  in  the  residual  amplitude  do  produce  any  noticeable  effects  in  the  bar  region. 


Figure  4-2:  The  ratio  R*  in  models  of  NGC  1073  where  the  residual  amplitude  Ar  = 
0 km^s“^kpc“^  (“best”  model,  filled  circles),  1000  km^s“^kpc“'  (open  circles),  and  2000 
km^s“^kpc“^  (plusses).  For  comparison,  the  primary  amplitude  A = 9000  km^s“^kpc“^ 
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Figure  4—2  shows  these  effects.  Here,  the  inclusion  of  any  residual  amplitude  inside 
the  bar  leads  to  worse  self-consistency. 

Variation  of  /q 

Figure  4-3  shows  the  results  of  varying  the  spiral  pitch  angle  /q.  Rather  large 
systematic  displacements  of  R*,  similar  to  those  due  to  variations  in  A (Figure  4-1), 
are  seen.  These  effects,  in  fact,  are  not  independent.  Contopoulos  and  Grosbpl  (1986) 
showed  that,  depending  on  the  assumed  geometry,  the  maximum  perturbed  density 
of  the  spiral  is  proportional  to  either  tan“^  io  (cylindrical  geometry)  or  Itan  io\~^  (flat 
geometry).  Decreasing  the  pitch  angle,  therefore,  produces  the  same  qualitative  effect  on 
R*  as  increasing  A.  The  plots  of  phase  difference  A9  exhibit  somewhat  larger  variations 
here  than  in  the  case  of  varying  A,  particularly  for  the  case  /’o  = -5°.  This  is  due  to 
the  fact  that,  given  the  same  perpendicular  offset  between  imposed  and  response  spiral, 
azimuthal  differences  between  the  two  are  magnified  more  in  the  case  of  a smaller 
pitch  angle. 

Variation  of  Es 
✓ 

A decrease  of  Cs  produces  a relative  increase  of  the  strength  of  the  spiral  outwards. 
As  can  be  seen  in  Figure  4-4,  this  causes  a greater  variation  in  the  imposed  spiral 
amplitude  than  in  the  response;  therefore,  R*  displays  a systematic  decrease  outwards 
when  the  value  of  Ss  is  too  small.  Conversely,  when  Ss  is  too  large,  the  imposed  spiral 
dies  off  too  quickly  and  the  ratio  R quickly  becomes  large.  These  results  hold  despite 
the  fact  that  the  value  of  A was  adjusted  in  order  to  keep  the  bar/spiral  perturbation 
amplitude  the  same  at  r = 2.5  kpc  in  all  three  cases.  The  phase  agreement  AO  is 
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approximately  the  same  for  the  two  lower  values  of  Es.  This  is  because  the  imposed 
spiral  remains  strong  enough  to  force  the  response  spiral  to  remain  almost  in  phase  with 
it.  When  Es  is  large,  though,  the  imposed  spiral  is  not  strong  enough  to  produce  good 
phase  agreement  beyond  corotation. 


Variation  of  ftp 


One  of  the  assumptions  of  our  modeling  technique  is  that  the  bar  figure  ends 
at  corotation.  This  assumption  is  based  on  results  of  orbit  calculations  in  realistic  bar 
potentials  (e.g.  Contopoulos  1980;  Teuben  and  Sanders  1985)  and  of  A^-body  simulations 
(e.g.  Sellwood  1981).  Still,  there  is  some  disagreement  as  to  the  exact  placement  of 
corotation,  and  many  authors  claim  that  bars  end  somewhat  inside  of  corotation.  In 
order  to  test  whether  either  of  these  changes  improve  self-consistency,  we  varied  the 
bar/spiral  pattern  speed  Qp  of  the  “best”  model  of  NGC  3992.  Figure  4—5  shows  the 
results  of  this  variation.  The  open  circles  represent  the  case  where  corotation  is  placed 
at  |a,  a being  the  bar  semimajor  axis.  The  filled  circles  represent,  as  usual,  the  results 
of  the  “best”  model.  The  plusses  represent  the  case  where  corotation  is  placed  at  |a. 
The  resonance  positions  marked  in  the  diagram  are  for  the  “best”  case  with  corotation 
at  r = a.  We  judge  that  the  best  results,  as  regards  both  R*  and  A9,  occur  when  the 
bar  is  assumed  to  end  at  corotation.  Still,  the  results  of  the  case  where  the  bar  ends 
before  corotation  are  not  very  much  worse  than  those  of  the  “best”  model.  Therefore, 
we  concur  with  those  who  claim  that  bars  end  at  or  slightly  inside  corotation,  and  in 
any  case  do  not  extend  beyond  corotation. 
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Figure  4-3:  The  ratio  R*  and  phase  difference  A9  in  models  of  NGC  3992  where  io  = -10° 
(“best”  model,  denoted  by  filled  circles),  io  = -5°  (open  circles),  and  io  = -15°  (plusses). 
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Figure  4^:  The  ratio  R*  and  phase  difference  A9  in  models  of  NGC  3992  where  Ss  = 0.2 
kpc~^  (open  circles,  A = 1213  km^s“^kpc“^),  0.4  kpc“^  (filled  circles,  “best”  model),  and  0.8 
kpc“^  (plusses,  A = 5437  km^s"^kpc"^). 
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Figure  4-5:  Amplitude  ratio  R*  and  phase  difference  in  models  of  NGC  3992  where  fip 
= 34.7  km  s“'kpc“^  (corotation  = 4a/3,  open  circles),  flp  = 43.6  km  s"^kpc~*  (“best”  model, 
corotation  = a,  filled  circles),  and  Dp  = 55.7  km  s“'kpc“*  (corotation  = 3a/4,  plusses). 
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Variation  of  A 


Figure  4—6  shows  the  results  of  varying  A,  the  z-thickness  of  the  bar/spiral 
perturbation.  We  see,  as  expected,  a similar  sort  of  effect  on  the  ratio  R as  varying 
the  amplitude  A directly  (Figure  4-1).  This  result  is  easily  understood.  Since  the  z- 
thickness  does  not  influence  the  potential,  or  therefore  the  orbital  behavior,  in  the  disk 
plane,  the  response  density  is  independent  of  A.  Decreasing  A,  then,  is  tantamount  to 
decreasing  the  amplitude  of  the  imposed  density  relative  to  the  response  (i.e.  increasing 
the  value  of  R*).  We  find  that,  for  all  three  galaxies  modeled,  the  optimum  value  of 
A is  in  the  range  from  0.05  kpc  to  0.2  kpc,  with  0.1  kpc  being  judged  as  yielding  the 
best  overall  self-consistency  in  each  case.  Finally,  since  A does  not  affect  the  response 
density,  the  phase  difference  AO  also  remains  unaffected. 


Figure  4-6:  The  ratio  R*  in  models  of  NGC  3992  where  A = 0.01  kpc  (open  circles),  0.1  kpc 
(filled  circles,  “best”  model),  and  1.0  kpc  (plusses). 
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Variation  of  62 

The  phase  angle  62  represents  the  angular  offset  between  the  innermost  part  of  the 
spiral  and  the  end  of  the  bar.  It  is  defined  such  that  positive  values  displace  the  spiral 
in  the  direction  of  pattern  rotation  (i.e.  the  spirals  “lead”  the  bar).  Negative  values 
correspond  to  spirals  that  “trail”  the  bar.  Figure  4—7  compares  the  results  of  setting 
62  = 15°  (open  circles)  and  30°  (plusses)  to  the  “best”  model  (filled  circles)  of  NGC 
3992.  We  note  large  deviations  in  both  R*  and  in  the  region  around  corotation, 
indicating  that  angle  offsets  between  the  ends  of  the  bar  and  the  starting  points  of  the 
spiral  arms  are  not  preferred  orientations  for  self-consistent  stellar  models  which  assume 
quasi-stationary  structure  that  rotates  at  a single  pattern  speed.  We  emphasize  positive 
values  of  O2  here  because  the  observations  of  our  sample  galaxies  seem  to  preclude  the 
existence  of  negative  offsets.  Nevertheless,  one  model  of  NGC  1398  was  calculated 
with  $2  = -15°.  The  results  of  that  run  are  as  unsuccessful  as  those  shown  in  Figure 
4-7,  but  the  deviations  of  R*  and  /S.6  are  in  the  opposite  sense.  We  had  hoped  for 
better  results  in  this  case  of  NGC  1073,  which  in  fact  displays  an  offset  of  some  30° 
between  the  bar  and  spirals  (Sandage  1961),  but  they  also  were  no  better  than  those 
shown  in  Figure  4—7  for  NGC  3992.  The  fact  that  NGC  1073  and  some  other  barred 
spirals  exhibit  angle  offsets  between  the  bars  and  spirals  may  indicate  either  that  the 
spirals  in  this  region  are  primarily  gaseous  in  nature  or  that  the  bar  and  spiral  strucures 
are  independent  entities  that  possess  separate  pattern  speeds.  There  is  some  numerical 
(e.g.  Sellwood  and  Sparke  1988)  and  observational  (e.g.  NGC  1365)  evidence  for 
independent  pattern  speeds  for  bars  and  spirals.  This  possibility  will  be  discussed  again 
in  Chapter  7. 
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Figure  4-7:  The  ratio  R*  and  phase  difference  A6  for  models  of  NGC  3992  in  which  62 
(“best”  model,  filled  circles),  15°  (open  circles),  and  30°  (plusses). 
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Variation  of  r\,  V2,  k\,  K2,  and  A4 


In  this  section  we  consider  the  effects  of  varying  the  remaining  parameters  of 
the  bar/spiral  perturbation.  These  are  the  cutoff  radii  (ri,r2)  and  steepnesses  (/ci,  K2) 
of  the  bar/spiral  and  the  amplitude  A4  of  the  A9  component.  Changes  in  all  of  these 
parameters  were  found  to  have  much  less  effect  on  model  self-consistency  than  those 
already  discussed. 

Figure  4—8  shows  the  effect  on  the  ratio  R*  of  the  “best”  model  of  NGC  3992  of 
varying  the  position  of  the  inner  bar/spiral  cutoff  radius  r\ . Not  surprisingly,  the  effects 
are  only  noticeable  in  the  vicinity  of  the  cutoff  itself  {r\  = 1.5  kpc  in  the  “best”  model  of 
NGC  3992).  Positioning  the  cutoff  closer  to  the  center  of  the  grid  (ri  = 0.5  kpc,  the  open 
circles  in  Figure  4-8)  reduces  somewhat  the  upward  divergence  of  R in  that  region. 


Figure  4-8;  The  ratio  R*  in  models  of  NGC  3992  where  r\  = 1.5  kpc  (“best”  model,  filled 
circles),  0.5  kpc  (open  circles),  and  2.5  kpc  (plusses). 
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In  fact,  since  the  amplitude  of  the  imposed  29  component  does  not  decrease  rapidly 
enough  toward  the  center  in  this  case,  R eventually  turns  downward  and  reaches  zero 
at  the  center.  This  case  is  arguably  an  improvement  over  the  “best”  model,  but  in  our 
opinion  the  overall  agreement  (i.e.  at  all  radii)  is  not  significantly  better. 

Figure  4—9  shows  the  results  of  varying  r\  in  the  case  of  NGC  1073.  Here  we  see 
much  a much  more  pronounced  effect  on  the  ratio  R than  in  the  case  of  NGC  3992  (or 
NGC  1398).  The  difference  is  that  the  cutoff  radius  is  relatively  farther  out  in  the  disk  in 
the  model  of  NGC  1073  than  it  is  in  the  models  of  either  of  the  other  two  galaxies.  Also, 
since  the  “best”  model  of  NGC  1073  employs  a rather  shallow  inner  cutoff  steepness 
Ki,  the  variation  in  r\  causes  a systematic  shift  in  R similar  to  the  effect  of  increasing 
the  bar/spiral  amplitude  A directly,  as  this  is  in  essence  what  has  happened.  In  the  case 
of  each  galaxy,  varying  r\  did  not  affect  the  phase  agreement  A9  at  all. 


Figure  4-9:  The  ratio  R*  in  models  of  NGC  1073  where  r\  = 2.95  (“best”  model,  filled  circles), 
2.0  kpc  (open  circles),  and  3.5  kpc  (plusses). 
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Models  of  NGC  3992  and  NGC  1073  were  also  calculated  in  which  the  position 
of  the  outer  cutoff  of  the  bar/spiral  perturbation,  r2,  was  varied  (an  outer  spiral  cutoff 
was  not  employed  in  the  case  of  NGC  1398).  These  models  showed  only  very  minor 
differences  with  the  “best”  models  of  each  case.  There  are  several  reasons  for  this.  Since 
V2  in  the  “best”  models  of  NGC  3992  and  NGC  1073  is  located  at  the  corresponding 
OLR,  the  spiral  amplitude  is  rather  small;  moreover,  the  cutoff  is  not  infinitely  sharp  but 
mediated  by  the  cutoff  steepness  k2  (1.0  and  0.5,  respectively,  for  the  “best”  models  of 
NGC  3992  and  NGC  1073,  denoting  rather  moderate  cutoff  steepnesses).  Therefore,  the 
exact  position  of  the  outer  cutoff  did  not  affect  the  self-consistencies  of  these  models 
to  a significant  extent. 

The  inner  and  outer  cutoff  steepnesses  ki  and  k2  were  varied  in  each  case  where 
the  corresponding  cutoff  was  employed.  Decreasing  either  value  (i.e.  making  the  cutoff 
too  shallow)  resulted  in  a “leakage”  of  excess  bar/spiral  perturbing  potential  beyond  the 
cutoff  and  a corresponding  depression  of  R*,  the  same  effect  seen  in  Figure  4—1  when 
A is  too  large.  The  opposite  effect  occurred  when  the  cutoff  steepnesses  were  taken  to 
be  too  large.  In  all  cases  the  phase  agreement  A0  was  essentially  unaffected. 

We  also  calculated  models  in  which  we  included  A9  components  with  amplitudes 
20%  and  40%,  respectively,  of  the  20  amplitude  A of  the  “best”  model.  We  find 
that  the  ratio  R remains  almost  completely  unaffected.  The  phase  agreement  A0 
is  altered  somewhat,  particularly  near  the  -4/1  resonance,  but  always  towards  larger 
negative  values  (i.e.  worse  self-consistency).  In  contrast,  Contopoulos  and  Grosbpl 
(1988)  found  that  the  inclusion  of  a moderate  AO  component  improves  self-consistency 
in  models  of  unbarred  spirals.  This  difference  may  be  due  in  part  to  the  relatively 
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large  fraction  of  stochastic  orbits  responsible  for  supporting  the  outer  spiral  structure 
in  the  barred  models. 

We  now  consider  variations  of  the  parameters  which  characterize  the  Ferrers  bar 
potential  and  the  exponential  disk.  In  all  cases  we  restrict  the  variations  to  within  the 
limits  set  by  the  surface  photometry  and  the  observed  rotation  curve,  and  we  adjust  the 
Plummer  sphere  halo  mass  M//  and  shape  parameter  bu  in  order  to  ensure  a good  fit 
to  the  observed  HI  rotation  curve. 

Variation  of  the  Bar  Semiaxes  a,  b,  and  c 

Since  varying  the  semimajor  axis  a of  the  Ferrers  bar  potential  is  essentially 
equivalent  to  varying  the  pattern  speed  Q,p  (the  main  effect  of  both  is  to  vary  the 
length  of  the  bar  relative  to  the  corotation  resonance),  it  is  not  surprising  that  the  effects 
on  R and  A9  in  each  case  are  qualitatively  the  same  as  those  shown  in  Figure  4—5. 
The  magnitude  of  the  changes  are  less,  though,  since  the  surface  photometry  constrains 
the  value  of  a (and  therefore  the  difference  between  a and  the  corotation  radius)  more 
tightly  than  was  allowed  when  ftp  was  varied.  Furthermore,  we  find  that  varying  b and 
c,  the  intermediate  and  minor  bar  axes,  within  the  limits  set  by  the  surface  photometry 
produces  essentially  no  effect  on  R*  and  A9. 

Variation  of  the  Bar  Mass  Mb 

Since  the  only  observational  constraints  on  the  masses  of  the  bars  of  our  program 
galaxies  are  the  (somewhat  uncertain)  inner  parts  of  the  corresponding  HI  rotation 
curves,  we  have  considerable  latitude  to  vary  Mg.  Figure  4—10  shows  the  effects  of 
varying  Mb  on  the  measures  of  self-consistency  in  the  case  of  NGC  3992.  We  see 
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that  R*  is  affected  the  most  in  the  region  of  the  bar,  while  Ad  is  affected  most  outside 
the  bar  in  the  region  of  the  spiral.  While  a moderate  decrease  in  bar  mass  improves 
the  self-consistency  as  regards  R*  in  the  region  of  the  bar  (e.g.  the  case  where  Mq 
= 7.5  X IO^Mq,  denoted  by  plusses  in  Figure  4—10),  this  improvement  is  more  than 
offset  by  the  deterioration  of  the  phase  agreement  in  the  outer  spiral.  We  consider  the 
model  with  the  highest  value  of  Mb  (i.e.  1.5  x 10^*^ M©)  to  be  the  most  self-consistent 
overall.  This  conclusion  is  strengthened  by  noting  that  Hunter  et  al.  (1988)  obtained  a 
comparable  mass  of  2.5  x 10^^ Mq  for  the  bar  of  NGC  3992  by  an  independent  modeling 
technique.  The  results  of  varying  Mb  in  the  cases  of  NGC  1073  and  NGC  1398  yielded 
results  similar  to  those  presented  here  for  NGC  3992. 


Variation  of  cq  and  eo 


Next  we  consider  the  central  surface  density  cq  and  inverse  scale  length  eo  of  the 
exponential  disk.  Figure  4—1 1 shows  the  effect  of  varying  cq  on  the  ratio  R in  models 
of  NGC  3992.  We  note  systematic  vertical  shifts  in  R , with  higher  values  of  central 
surface  density  corresponding  to  higher  values  of  response-to-imposed  amplitude  ratio. 
This  is  due  to  the  fact  that  when  cq  is  increased  and  the  amplitude  of  the  imposed  bar 
and  spiral  perturbation  are  held  fixed,  the  amplitude  of  the  bar/spiral  relative  to  the 
axisymmetric  background  is  decreased.  The  response  amplitude  does  not  suffer  such 
a decrease  because  the  amplitude  of  the  response  bar  and  spiral  is  modulated  by  cq 
through  the  weighting  of  the  orbits.  The  ratio  R*,  then,  tends  to  go  up  when  cq  is 


increased  and  down  when  it  is  decreased. 
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Figure  4-10:  The  ratio  R*  and  phase  difference  AO  for  models  of  NGC  3992  in  which  Mb  = 
1.5  X 1O^°M0  (“best”  model,  filled  circles),  1.5  x IO^Mq  (open  circles),  and  7.5  X 10^ M0 
(plusses). 
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Figure  4-11:  The  ratio  R*  in  models  of  NGC  3992  where  cq  = 75OM0/pc^  (“best”  model, 
filled  circles),  SOOMq/pc^  (open  circles),  and  lOOOM0/pc^  (plusses). 

Figure  4-12  shows  the  effect  of  varying  the  inverse  disk  scale  length  eo  on  the 
ratio  R*  in  models  of  NGC  3992.  Here  again  we  see  systematic  shifts  of  the  response- 
to-imposed  amplitude  ratio  and  for  the  same  reasons  outlined  above  for  the  variation 
of  central  surface  density  cq.  A lower  value  of  eo  means  a less  rapid  dropoff  in 
axisymmetric  density  away  from  the  center.  This,  in  turn,  results  in  a progressive 
decrease  with  radius  of  the  bar/spiral  amplitude  relative  to  underlying  disk,  and  a 
corresponding  increase  in  R*.  The  opposite  is  true  for  the  case  of  higher  sq.  We 
can  clearly  see  this  divergence  in  R away  from  center  for  the  different  cases. 

We  detect  essentially  no  changes  in  the  phase  agreement  A0  when  cq  and  £q 
are  varied.  This  result  might  seem  surprising  at  first,  but  it  really  is  not.  The  phase 
agreement  between  the  response  and  imposed  bar/spiral  is  basically  determined  by  the 
shapes  of  the  orbits,  which,  in  turn,  are  determined  by  the  total  potential  in  the  disk. 
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We  noted  above  that  we  altered  the  halo  parameters  at  the  same  time  as  we  changed  the 
exponential  disk  parameters  in  order  to  maintain  a good  fit  to  the  observed  rotation  curve. 
Therefore,  we  have  compensated  the  gravitational  potential  in  the  disk  for  changes  in 
the  disk  parameters  and  ensured  that  the  orbits  have  maintained  their  essential  shapes. 
This  is  why  the  is  not  affected  by  variations  in  the  disk  parameters. 


Figure  4-12:  The  ratio  R*  in  models  of  NGC  3992  where  eo  = 0.235  kpc  ^ (“best”  model, 
filled  circles),  0.157  kpc“^  (open  circles),  and  0.353  kpc“*  (plusses). 

Variation  of  cto  and  <7r 


Finally,  we  consider  the  effects  on  model  self-consistency  of  varying  the  central 
value  (Jo  and  radial  slope  Cr  of  the  velocity  dispersion  employed.  Following  the  results 
of  observations  (e.g.  Jarvis  et  al.  1988)  and  of  V-body  simulations  (e.g.  Sparke  and 
Sellwood  1987;  Pfenniger  and  Friedli  1991),  we  consider  here  only  those  models  in 
which  the  velocity  dispersion  decreases  with  radius.  We  have  calculated  models  in 
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which  the  velocity  dispersion  is  constant  across  the  disk,  but  these,  in  general,  yield 
results  inferior  to  those  in  which  the  dispersion  decreases.  Figure  4—13  shows  the 
effects  of  varying  ao  and  ar  on  R*  and  A9  in  models  of  NGC  3992.  We  see  that 
a lower  velocity  dispersion  (open  circles  in  Fig.  4—13)  leads  to  worse  agreement  as 
regards  R*  but  good  agreement  as  regards  A9.  On  the  other  hand,  a higher  velocity 
dispersion  (plusses  in  Fig.  4—13)  yields  slightly  worse  agreement  in  both  R*  and  A^. 
These  results  are  understood  by  realizing  that  lower  velocity  dispersions  lead  to  a more 
sharply  peaked  response  (and  hence  a greater  29  response  component)  which  lies  close 
to  the  imposed  potential  minima.  Higher  velocity  dispersions,  conversely,  tend  to  smear 
out  the  response  (and  hence  lower  the  29  response  component),  as  well  as  cause  it  to 
depart  farther  from  the  minima  of  the  potential.  Thus  it  seems  that  the  appropriate 
velocity  dispersion  is  a compromise  between  these  two  competing  effects.  The  velocity 
dispersion  which  works  best  for  our  model  of  NGC  3992  has  a central  value  of  100  km 
s“^  and  drops  to  around  25  km  s“^  near  the  OLR  (10.6  kpc). 

Interpretation  of  Results 

In  this  final  section  we  briefly  summarize  the  conclusions  that  can  be  drawn  from 
the  results  of  varying  the  model  parameters.  The  properties  which  seem  to  exert  the 
most  influence  on  self-consistency  are  (1)  the  bar/spiral  amplitude  A,  (2)  the  spiral  pitch 
angle  io,  (3)  the  inverse  spiral  scale  length  Ss,  (4)  the  perturbation  thickness  A,  (5)  the 
bar/spiral  pattern  speed  (6)  the  bar  mass  Mb,  (7)  the  cental  surface  density  Cq  and 
inverse  scale  length  eo  of  the  exponential  disk,  and  (8)  the  central  value  ao  and  radial 
slope  ar  of  the  velocity  dispersion.  Parameters  which  control  the  detailed  structure  of 
the  bar/spiral  amplitude  such  as  the  cutoff  radii  r\  and  V2  and  steepnesses  k\  and  K2,  as 
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kpc 

Figure  4-13:  The  ratio  R*  and  phase  difference  A0  in  models  where  ao  = 100  km  s“^  = 

-7.0  km  s“^kpc“^  (“best”  model,  filled  circles),  ao  = 30  km  s“^  a^  = -1.0  km  s”^kpc"^  (open 
circles),  and  ao  = 150  km  s“^  a^  = -11.0  km  s"^kpc“^  (plusses). 
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well  as  the  residual  amplitude  Ar,  do  exert  influence  on  model  self-consistency,  but  their 
effects  are  secondary  to  the  parameters  which  control  its  global  shape  and  strength  (i.e. 
A,  io,  Ss,  A,  and  Mb).  Also,  while  the  bar  length  a affects  self-consistency  significantly, 
this  simply  reflects  the  influence  of  the  pattern  speed  flp.  Thus  we  may  group  the 
important  parameters  into  four  main  categories:  (i)  those  that  control  the  global  shape 
and  amplitude  of  the  bar/spiral  potential,  (ii)  the  pattern  speed,  (iii)  those  that  control 
the  magnitude  and  shape  of  the  azimuthally  averaged  radial  mass  distribution,  and  (iv) 
those  that  control  the  magnitude  and  shape  of  the  velocity  dispersion.  Finally,  we  again 
note  the  rather  strong  preference  in  the  models  of  this  thesis  for  the  spiral  arms  to 
emanate  directly  from  the  ends  of  the  bar. 

Since  we  only  have  three  galaxies  in  our  sample,  it  is  difficult  to  discern  possible 
correlations  of  the  parameters  of  the  “best”  models  with  Hubble  type.  Still,  some 
comments  are  warranted.  Three  basic  classification  criteria  for  spiral  (and  therefore 
barred  spiral)  galaxies  are  given  by  Sandage  (1961):  (1)  the  openness  of  the  spiral 
arms,  (2)  the  degree  of  resolution  of  the  arms  into  stars,  and  (3)  the  relative  size  of 
the  unresolved  nuclear  region.  Although  there  are,  of  course,  exceptions,  earlier  type 
galaxies  (here  early,  intermediate,  and  late  type  barred  spirals  refer  to  SBa,  SBb,  and 
SBc,  respectively,  in  Hubble’s  [1926]  original  classification  scheme)  tend  to  have  both 
smoother,  more  tightly  wound  spiral  arms  that  are  less  resolvable  into  individual  stars 
and  relatively  larger  unresolved  nuclear  regions  (i.e.  “bulges”).  The  distribution  of 
spiral  pitch  angles  in  our  “best”  models  is  consistent  with  these  expectations.  The 
“best”  model  of  our  representative  early-type  galaxy  NGC  1398  has  a pitch  angle  of 
-6°,  while  the  “best”  models  of  the  two  later-type  representatives  (NGC  3992  and  NGC 
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1073)  both  have  i‘o  = -10°.  While  we  might  have  expected  the  “best”  model  of  NGC 
1073  to  display  a somewhat  larger  pitch  angle,  we  note  that  the  observed  scatter  in 
pitch  angle  as  a function  of  Hubble  type  is  rather  large  (Kennicutt  1981)  and  that  the 
value  of  i'o  for  NGC  1073  is  not  unusual.  We  are  unable  to  say  anything  based  on  our 
model  results  concerning  the  latter  two  criteria.  This  is  because  criterion  (2)  is  uniquely 
observational  in  nature  and  criterion  (3)  refers  to  components  (bulges)  not  included  in 
the  present  models. 

One  final  point  we  would  like  to  make  concerns  an  observation  of  Elmegreen  and 
Elmegreen  (1985).  In  their  paper  they  distinguish  between  two  types  of  bars,  and 
exponential,  depending  on  the  shape  of  the  luminosity  profile  along  the  major  axis. 
They  find  that  early-type  barred  spirals  are  more  likely  by  far  to  possess  a flat  bar, 
while  exponential  bars  are  preferentially  found  in  late-type  galaxies.  The  results  of  our 
modeling  tend  to  corroborate  this  finding.  In  particular,  the  “best”  models  of  NGC 
3992  and  NGC  1398  require  bar/spiral  potentials  which  extend  significantly  inwards 
of  corotation,  while  the  inner  cutoff  radius  of  the  bar/spiral  in  the  “best”  model  of 
NGC  1073  is  located  at  corotation.  At  the  same  time  we  note  that  our  assumed  Ferrers 
bar  potential  more  closely  represents  an  exponential  bar  than  a flat  bar.  We  interpret 
the  inward  extensions  of  the  bar/spiral  perturbations  in  the  models  of  NGC  3992  and 
NGC  1398  as  compensating  for  the  intrinsic  “non-flatness”  of  the  Ferrers  bar  potential. 
Combes  and  Elmegreen  (1992)  give  further  results  concerning  flat  versus  exponential 
bars.  Their  A-body  simulations  indicate  that  the  flat  bars  of  early-type  galaxies  are 
limited  in  extent  by  their  corotation  radii  (i.e.  by  stochasticity),  whereas  the  exponential 
bars  of  late-type  galaxies  are  limited  by  the  scale  length  of  the  disk  in  which  they  are 
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located  (i.e.  the  bar  stops  growing  when  the  matter  density  in  between  corotation  and 
the  OLR  is  not  sufficient  to  absorb  the  angular  momentum  transported  outwards  by  the 
growing  bar).  They  claim,  therefore,  that  the  ends  of  exponential  bars  do  not  necessarily 
correspond  to  any  major  resonance  radius.  We  did  not  check  this  possibility  here. 


CHAPTER  5 
STOCHASTIC  ORBITS 


In  this  chapter  we  examine  the  role  of  stochasticity  in  the  models  presented  in 
Chapter  3.  Again,  since  we  are  considering  here  only  the  surface  density  response, 
we  restrict  our  attention  to  stochastic  motion  in  the  disk  plane.  While  this  restriction 
simplifies  the  present  analysis,  it  completely  misses  potentially  important  phenomena 
such  as  vertical  resonances,  Arnold  diffusion,  complex  instability,  and  collisions  of 
bifurcations,  all  of  which  only  appear  in  cases  of  at  least  three  degrees  of  freedom 
(Contopoulos  1987).  Indeed,  analyses  of  the  orbital  structures  of  realistic  bar  potentials 
(e.g.  Pfenniger  1984a,  1985;  Pfenniger  and  Friedli  1991)  and  the  results  of  3-D  A-body 
simulations  (e.g.  Combes  and  Sanders  1981;  Pfenniger  and  Friedli  1991)  indicate  that 
vertical  instability  strips  in  the  periodic  orbit  families,  as  well  as  global  buckling  (“fire 
hose”)  instabilities  of  the  bars  themselves,  can  play  significant  roles  in  determining 
the  final  3-D  (often  box  or  peanut)  shape  that  a given  bar  assumes.  These  vertically 
thickened  bars  of  the  simulations,  in  fact,  have  been  identified  with  the  similarly  shaped 
bulges  seen  in  some  edge-on  galaxies  (Combes  and  Sanders  1981). 

These  results  notwithstanding,  we  are,  for  several  reasons,  of  the  opinion  that  the 
2-D  problem  should  be  addressed  in  more  detail  first.  For  one,  the  collective  long- 
term surface  density  response  of  stochastic  orbits  in  realistic  two-dimensional  barred 
spiral  potentials  has  simply  never  been  tested.  Pfenniger  (1984b),  in  his  attempts  to 
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construct  self-consistent  2-D  bars,  found  that  he  usually  had  to  include  roughly  10% 
stochastic  orbits  in  order  to  achieve  success.  Other  authors  (e.g.  Teuben  and  Sanders 
1985)  make  the  claim  that  if  an  orbit  is  stochastic  (i.e.  it  only  respects  one  integral  of 
motion,  the  energy),  then  its  time-averaged  density  response  fills  more  or  less  uniformly 
the  region  bounded  by  its  zero-velocity  curve.  However,  Contopoulos  and  Polymilis 
(1993)  have  shown  that  in  near-integrable  systems  (such  as  galactic  models),  stochastic 
orbits  may  remain  for  long  times  in  the  phase-space  neighborhood  of  unstable  periodic 
orbits.  Contopoulos  and  Kaufmann  (1992)  have  emphasized  this  “stickiness”  property 
as  the  phenomenon  that  explains  both  the  almost-regular  behavior  of  stochastic  orbits 
over  long  times  and  the  fact  that  these  orbits  do  not  escape  to  infinity  for  long  times, 
even  though  their  energy  exceeds  the  escape  energy  and  their  motion  is  not  restricted 
by  any  other  integral  of  motion. 

A second  reason  to  study  the  density  response  of  stochastic  orbits  in  realistic  two- 
dimensional  galactic  potentials  is  indicated  by  the  results  of  Sparke  and  Sellwood’s 
(1987)  dissection  of  a 2-D  A-body  bar.  First,  they  found  that  the  bar  was  comprised 
of  orbits  trapped  and  semi-trapped  around  the  central  x\  family  of  periodic  orbits. 
The  trapped  orbits  are  true  quasi-periodic  orbits  which  lie  on  invariant  tori  in  phase 
space.  The  semi-trapped  orbits  are  actually  stochastic  orbits  whose  exploration  of  the 
available  phase  space  is  restricted  for  long  times  by  cantori  (invariant  tori  which  have 
developed  an  infinity  of  “holes”  due  to  an  increase  in  the  perturbation,  thereby  allowing 
communication  between  formerly  well-separated  regions  of  stochasticity).  We  discuss 
these  ideas  in  more  detail  below.  Secondly,  and  perhaps  more  importantly  for  barred 
spiral  models,  they  found  that  there  existed  a substantial  “hot”  population  of  particles 
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which  had  sufficient  energy  to  traverse  the  bar  and  the  outer  disk  (i.e.  their  Jacobi 
constants  exceeded  the  value  at  the  Lagrange  points  L\  and  L2).  Pfenniger  and  Friedli 
(1991)  also  found  such  a “hot”  population  in  their  similar  study  of  a 3-D  A^-body  bar. 
As  we  shall  see,  it  is  this  class  of  orbits  which  plays  the  crucial  role  in  supporting  the 
imposed  spiral  structure  of  our  models  from  the  ends  of  the  bar  (near  corotation)  to  a 
point  in  the  outer  disk  where  regular  orbits  again  can  dominate  the  spiral  structure. 

Finally,  observations  of  edge-on  spirals  seem  to  indicate  that  at  least  some  bars 
are  as  thin  as  the  rest  of  the  disk.  This  statistical  argument  says  that  if  bars  are  as 
common  in  edge-on  systems  as  in  all  disk  galaxies  (-25  - 35%),  then  we  should  see  a 
substantial  fraction  of  thick  disks  if  bars  are  substantially  thicker  than  the  disks  in  which 
they  reside.  Because  we  do  not,  bars  must  be  as  thin  as  the  disks  of  their  host  galaxies. 
While  the  observation  of  box  and  peanut-shaped  bulges  tends  to  undercut  this  argument, 
the  proportion  of  box  and  peanut-shaped  bulges  in  edge-on  galaxies  (-20%,  e.g.  Shaw 
1987)  does  not  match  the  proportion  of  bars  in  all  disk  galaxies.  Also,  Athanassoula 
(cf.  Sell  wood  and  Wilkinson  1992)  raises  the  concern  that  the  sizes  of  the  observed 
box/peanut-shaped  bulges  relative  to  the  extent  of  the  disk  are  less  than  the  sizes  of 
strong  bars  to  their  host  galaxies.  These  facts  imply  that  a significant  proportion  of 
bars  are  thin  as  well.  If  this  is  true,  the  importance  of  phenomena  which  depend  on  a 
strong  coupling  between  motion  in  the  disk  and  motion  in  the  third  dimension  will  be 
lessened.  A typical  value  quoted  for  the  ratio  of  bar  thickness  to  bar  major  axis  is  1:10 
(e.g.  Kormendy  1982;  Wakamatsu  and  Hamabe  1984). 

A number  of  studies  have  addressed  the  problem  of  stochasticity  in  two-dimensional 
galactic  models  (e.g.  Contopoulos  1981,  1983,  1987;  Contopoulos  and  Grosbpl  1989; 
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Athanassoula  et  al.  1983;  Teuben  and  Sanders  1985).  The  transition  from  the  com- 
pletely ordered  motion  of  an  integrable  system  (i.e.  a system  in  which  all  particles  obey 
the  same  number  of  integrals  of  motion  as  there  are  degrees  of  freedom)  to  the  stochastic 
motion  of  a chaotic  system  is  now  fairly  well-understood.  The  introduction  to  an  inte- 
grable two-dimensional  system  of  a small  perturbation  which  couples  the  independent 
motions  renders  the  system  near-integrable.  That  is,  the  system  contains  an  infinite 
number  of  stable  and  unstable  periodic  orbits,  and  regions  of  chaos  are  found  in  the 
vicinity  of  each  unstable  orbit.  Because  the  perturbation  is  small,  however,  the  regions 
of  chaos  are  also  small,  and  the  system  retains  a large  set  of  ordered  motions.  These 
motions  lie,  in  phase  space,  on  toroidal  surfaces  called  KAM  (Kolmogorov,  Arnold  and 
Moser,  e.g.  Moser  1973)  tori,  which  separate  the  chaotic  regions.  The  intersections  of 
the  KAM  tori  by  a Poincare  surface  of  section  (cf.  Poincare  1892)  are  closed  invariant 
curves.  As  the  perturbation  increases,  the  relative  sizes  of  the  stochastic  regions  also 
increase.  There  is,  however,  no  communication  between  different  stochastic  regions 
until  the  perturbation  exceeds  a critical  value  and  the  so-called  last  KAM  curve  sep- 
arating them  is  destroyed.  The  evolution  of  the  last  KAM  curve  has  been  described 
by  Greene  (1979)  and  Shenker  and  Kadanoff  (1982).  When  the  last  KAM  curve  is 
destroyed,  it  develops  an  infinity  of  holes  which  form  a Cantor  set;  therefore,  it  is 
called  a cantorus  (Aubry  1978;  Percival  1979;  Mather  1982;  Katok  1982;  Aubry  and 
Le  Doeron  1983).  The  existence  of  these  holes  allows  diffusion  of  orbits  through  the 
cantorus  and  establishes  communication  between  previously  well-separated  regions  of 
chaos.  The  rate  of  diffusion  through  the  holes  of  the  cantorus  for  perturbations  slightly 
above  the  critical  value  has  been  estimated  theoretically  by  Bensimon  and  Kadanoff 
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(1984)  and  by  MacKay,  Meiss,  and  Percival  (1984).  This  rate  depends  very  much  on 
the  sizes  of  the  holes  in  the  cantorus,  and  therefore  on  the  amplitude  of  the  perturbation 
(Contopoulos  and  Barbanis  1989).  This  explains  why  stochastic  orbits  often  do  not 
cover  the  available  phase  space  in  a smooth  way.  Also,  if  the  diffusion  time  scale  is 
much  greater  than  the  average  orbit  period  (as  may  be  the  case  in  galactic  models), 
then  we  may  consider,  for  practical  purposes,  a cantorus  to  be  an  effective  barrier  to 
the  diffusion  of  stochastic  orbits. 

In  the  next  section  we  present  several  Poincare  surfaces  of  section  (each  repre- 
senting a different  value  of  the  Jacobi  constant)  for  the  most  successful  models  of  NGC 
3992,  NGC  1073,  and  NGC  1398  that  were  presented  in  Chapter  3.  The  surface  of 
section  is  a very  useful  tool  not  only  for  visualizing  the  periodic  orbits  that  are  present 
at  a given  value  of  the  Jacobi  constant,  but  also  for  gauging  the  amount  of  trapping 
done  by  each  periodic  orbit  and  the  amount  of  stochasticity  present.  We  construct  the 
surface  of  section  by  numerically  integrating  an  orbit  and  plotting  one  position  and  and 
one  velocity  component  whenever  it  crosses  some  reference  axis  in  a specified  direc- 
tion. For  example,  in  our  surfaces  of  section  we  take  the  reference  axis  to  be  along  the 
intermediate  axis  of  the  bar  (i.e.  y = 0),  and  plot  x and  x whenever  the  orbit  crosses 
the  reference  axis  with  y > 0.  Such  plotted  points  are  called  consequents.  We  shall 
see  how  the  orbital  structure  of  each  model,  with  emphasis  on  the  stochastic  regions, 
evolves  as  a function  of  the  Jacobi  constant.  After  presenting  the  surfaces  of  section, 
we  shall  examine  briefly  the  qualitative  behavior  of  the  stochastic  orbits  as  a function  of 
Jacobi  constant.  We  conclude  the  chapter  with  estimates  of  the  percentage  of  stochastic 


orbits  in  each  “best”  model. 
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Surfaces  of  Section 

In  Figures  5-1  through  5-9  we  present  surfaces  of  section  showing  the  evolution  of 
orbital  structure  and  the  onset  of  stochasticity  in  the  models  of  Chapter  3.  Each  surface 
of  section  is  generated  by  plotting  100  consequents  per  orbit  for  a large  number  of 
orbits  (about  100  per  surface  of  section),  all  with  the  same  value  of  the  Jacobi  constant. 
Figures  5-1  through  5-3  refer  to  the  model  of  NGC  3992,  Figures  5-4  through  5-6 
the  model  of  NGC  1073,  and  Figures  5-7  through  5-9  the  model  of  NGC  1398.  The 
limiting  curves  that  delineate  accessible  regions  are  shown  as  dashed  lines.  For  a given 
model  the  sequence  of  figures  is  arranged  in  order  of  increasing  Jacobi  constant. 

Figure  5-1  shows  surfaces  of  section  for  the  two  lowest  values  of  the  Jacobi 
constant  Ej  considered  in  the  model  of  NGC  3992.  The  most  striking  feature  in  the  top 
panel  of  Figure  5-1  {Ej  = -213438  km^  s“^)  is  the  predominant  regularity  of  the  orbits, 
both  in  the  region  of  the  bar  and  in  the  outer  disk.  (The  consequents  of  orbits  in  the 
bar  are  found  in  the  bounded  region  centered  on  x = 0,  while  those  in  the  outer  disk  are 
found  in  the  two  accessible  regions  on  the  left  and  right  extremes  of  the  diagram.  Also, 
since  we  focus  on  the  behavior  of  orbits  in  our  models,  we  include  only  consequents 
generated  by  direct  orbits.  That  is  why  the  consequents  in  the  bar  region  generally  have 
positive  X values  while  those  in  the  outer  disk  generally  have  negative  x values.  We  have 
included,  in  the  bar  region  of  the  top  panel  of  Figure  5-1,  one  invariant  curve  centered 
near  [x  = -1  kpc,  i:  = 0 km  s~^]  generated  by  a retrograde  orbit  to  show  that  there  is 
indeed  structure  to  these  regions.)  At  this  lowest  value  of  Ej,  which  corresponds  in  the 
axisymmetric  case  to  a circular  radius  Kc  = 12.0  kpc  (also  Cc  = 1.68  kpc),  the  orbital 
structure  in  the  bar  region  and  in  the  outer  disk  region  is  dominated  by  the  x\  family. 
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Figure  5-1:  Surfaces  of  section  in  the  model  of  NGC  3992.  Top:  Jacobi  constant 
equals  -213438  km^  s”^.  Bottom:  Jacobi  constant  equals  -202713  km^  s“^. 
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Figure  5-2:  More  surfaces  of  section  in  the  model  of  NGC  3992.  Top:  Jacobi 
constant  equals  -198084  km^  s"^.  Bottom:  Jacobi  constant  equals  -194202  km^  s~^. 
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Figure  5-3:  More  surfaces  of  section  in  the  model  of  NGC  3992.  Top:  Jacobi 
constant  equals  -191307  km^  s~^.  Bottom:  Jacobi  constant  equals  -189692  km^  s~^. 
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Figure  5^:  Surfaces  of  section  in  the  model  of  NGC  1073.  Top:  Jacobi  constant 
equals  -34143  km^  s“^.  Bottom:  Jacobi  constant  equals  -32393  km^  s“^. 
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Figure  5-5:  More  surfaces  of  section  in  the  model  of  NGC  1073.  Top:  Jacobi 
constant  equals  -30871  km^  s“^.  Bottom:  Jacobi  constant  equals  -29746  km^  s~^. 
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Figure  5-6:  More  urfaces  of  section  in  the  model  of  NGC  1073.  Top\  Jacobi 
constant  equals  -29405  km^  s“^.  Bottom:  Jacobi  constant  equals  -29265  km^  s“^. 
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Figure  5-7:  Surfaces  of  section  in  the  model  of  NGC  1398.  Top:  Jacobi  constant 
equals  -314204  km^  s“^.  Bottom:  Jacobi  constant  equals  -306600  km^  s“^. 
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Figure  5-8:  More  surfaces  of  section  in  the  model  of  NGC  1398.  Top:  Jacobi 
constant  equals  -299940  km^  s“^.  Bottom:  Jacobi  constant  equals  -294526  km^  s“^. 
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Figure  5-9:  More  surfaces  of  section  in  the  model  of  NGC  1398.  Top:  Jacobi 
constant  equals  -290724  km^  s“^.  Bottom:  Jacobi  constant  equals  -288987  km^  s“^. 
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The  xi  orbits  at  this  value  of  £7  are  highly  elongated  along  the  bar  in  the  bar  region 
(see  Figure  3-6)  and  essentially  circular  in  the  outer  disk.  Since  this  outer  branch  of  the 
Xi  family  is  located  beyond  the  OLR  where  the  matter  density  is  low,  it  is  of  negligible 
importance  in  our  models.  The  lower  panel  of  Figure  5-1  represents  the  case  where 
Ej  = -202713  km^  s“^  (rc  = 10  kpc,  2.44  kpc).  The  orbital  structure  of  the  bar  is  still 
dominated  by  the  Xi  family,  while  outside  we  see  two  main  stable  orbit  families  and 
a significant  amount  of  stochasticity.  The  two  stable  families  represent  -2/1  resonant 
orbits  near  the  OLR.  The  family  with  deformed  and  elongated  invariant  curves  centered 
near  (x  = -9.2  kpc,  ;fc  = 0 km  s“^)  is  the  same  orbit  family  as  the  X\  family  in  the  top 
panel  of  Figure  5-1.  Instead  of  being  circular,  though,  at  this  value  of  Ej  the  orbits 

0 

of  this  family  are  elongated  in  the  y direction  and  are  out  of  phase  with  the  imposed 
spiral.  Orbits  of  the  other  -2/1  family,  centered  near  (x  = -11  kpc,  jc  = 10  km  s~^), 
are  elongated  along  the  x direction  and  are  in  phase  with  the  imposed  spiral  (note:  the 
consequents  outside  the  bar  with  x > 0 kpc,  seen  in  the  lower  panel  of  Figure  5-1  and 
later  plots,  are  formed  by  loops  counter  to  the  general  motion  of  the  outer  orbits  around 
the  galaxy.  Figure  5-10  shows  schematically  how  such  consequents  are  generated.). 

Figure  5-2  shows  surfaces  of  section  for  Ej  - -198084  km^  s“^  {Tc  = 9 kpc,  2.94 
kpc;  top  panel)  and  Ej  = -194202  km^  s“^  {Tc  = 8 kpc,  3.52  kpc;  bottom  panel).  At  the 
lower  value  of  Ej  we  continue  to  see  almost  complete  regularity  as  regards  the  orbits 
in  the  bar.  Outside  the  bar,  except  for  a small  stable  -2/1  family  near  (x  = -9.2  kpc,  x 
= 0 km  s“^),  stochasticity  is  dominant.  While  there  are  stable  regions  evident  at  large 
values  of  x and  Ixl,  their  dynamical  significance  is  small  because  the  amount  of  matter 
existing  in  these  regions  of  phase  space  is  small  (the  matter  density  in  these  surfaces  of 
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section  decreases  exponentially  with  increasing  x and  Lxl).  In  the  lower  panel  of  Figure 
5-2  we  see  a significant  change  in  the  orbital  structure  of  the  bar.  At  this  value  of  Ej 
the  xi  family  in  the  bar  has  been  broken  by  a gap  into  two  4/1  resonant  families,  located 
near  (x  = 0.2  kpc,  A:  = 0 km  s“^)  and  (x  = 2.2  kpc,  x = 0 km  s~^).  The  first  4/1  family, 
however,  is  dynamically  insignificant  because  it  is  stable  over  only  a very  small  range 
of  energy.  Outside  the  bar  region  there  are  two  regions  of  stability.  The  one  centered 
roughly  at  (x  = -9  kpc,  j:  = 10  km  s“^)  is  an  extension  of  the  -2/1  family  of  the  top 
panel  of  Figure  5-2.  The  other,  located  approximately  at  (x  = -7.8  kpc,  x = -20  km 
s“^),  is  a -4/1  resonant  family  that  is  stable  only  over  a small  range  of  values. 


Figure  5-10:  Schematic  drawing  showing  how  loops  in  the  orbits  outside  of  corotation  generate 
consequents  on  the  surface  of  section.  A portion  of  an  orbit  containing  such  a loop  is  shown, 
together  with  arrows  indicating  the  direction  of  motion  along  the  orbit.  The  drawn  circle 
represents  corotation,  and  the  directions  of  the  general  motion  are  also  indicated  with  arrows. 

Finally,  in  Figure  5-3  we  see  surfaces  of  section  for  E/  = -191307  km^  s~^  (rc  = 

1 kpc,  4.21  kpc;  top  panel)  and  Ej  = -189692  km^  s“^  {Vc  = 6 kpc,  5.04  kpc;  bottom 
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panel).  These  two  surfaces  of  section  are  extremely  similar.  Inside  the  bar  region 
the  stochasticity  is  almost  complete.  Outside  the  bar  there  is  also  a large  region  of 
stochasticity.  Regions  of  regular  motion  are  confined  to  large  x values  (relative  to  r^) 
and  large  iJcl  values  where  matter  density  is  low,  thus  their  dynamical  significance  in 
this  range  of  Ej  is  small.  In  the  bottom  panel  of  Figure  5-3  the  consequents  around  (x 
= 3.8  kpc,  X = 0 km  s“^)  are  generated  by  stochastic  orbits  that  traverse  the  Lagrange 
points  Li  and  L2  and  travel  between  the  bar  and  the  outer  disk.  These  orbits,  which 
only  exist  for  values  of  Ej  > -191188  km^  s“^  (i.e.  the  Jacobi  constant  at  Li  and  L2, 
form  the  so-called  “hot”  population.  As  we  shall  see  shortly,  they  are  instrumental  in 
supporting  the  imposed  spiral  in  the  region  near  corotation. 

Figures  5-4  through  5-6  show  surfaces  of  section  derived  from  the  model  of  NGC 
1073  that  was  presented  in  Chapter  3.  Figures  5-7  through  5-9  show  analogous  figures 
for  the  case  of  NGC  1398.  In  both  progessions  of  plots  we  see  an  evolution  similar  to 
that  in  the  case  of  NGC  3992.  At  relatively  low  values  of  Ej  the  orbital  structure  within 
the  bar  is  very  regular.  As  the  Jacobi  constant  is  increased,  however,  a large  degree  of 
stochasticity  develops.  The  stochasticity  increases  with  increasing  Ej,  and  it  dominates 
the  motion  in  the  bar  region  when  Ej  is  sufficiently  close  to  the  Jacobi  constant  at  L\ 
and  L2.  Stochasticity  in  the  bar  region  of  the  model  of  NGC  1073  develops  somewhat 
more  rapidly  and  completely  than  in  the  model  of  NGC  1398.  This  is  due  primarily 
to  the  high  value  of  the  bar  axial  ratio  alb  in  the  model  of  NGC  1073.  Another  point 
of  difference  is  the  lack  of  any  significant  regular  motion  in  the  outer  disk  of  NGC 
1398  at  low  values  of  Ej.  This  is  due  to  the  fact  that,  at  large  distances,  the  outer 
spiral  in  the  model  of  NGC  1398  is  relatively  stronger  than  in  the  model  of  NGC  1073, 
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and  it  induces  greater  stochasticity  at  the  lower  values  of  Ej.  The  large  regions  of 
stochasticity  and  the  few  significant  periodic  orbit  families  that  exist  in  these  models, 
and  in  the  model  of  NGC  3992  as  well,  indicate  why  we  chose  to  include  the  “circular” 
orbits  mentioned  in  Chapters  2 and  3.  Since  we  needed  some  means  of  estimating  the 
density  response  of  the  orbits  in  these  predominantly  stochastic  regions,  and  because 
this  density  response  is  not  very  sensitive  to  the  particular  choice  of  initial  conditions, 
we  chose  to  start  the  orbits  with  the  initial  conditions  of  the  truly  circular  orbits  of 
the  axisymmetrized  case.  In  the  next  section  we  investigate  the  various  behaviors  of 
stochastic  orbits  in  the  configuration  space  of  the  disk  plane. 

Individual  Stochastic  Orbits 

In  this  section  we  examine  the  various  types  of  generic  stochastic  behavior 
exhibited  by  orbits  in  our  models.  Since  these  types  of  behavior  are  the  same  in 
all  three  models,  we  shall  restrict  our  attention  to  illustrative  examples  drawn  from  the 
model  of  NGC  3992  that  was  presented  in  Chapter  3.  Basically,  there  are  four  distinct 
types  of  behavior  that  a stochastic  orbit  can  display.  The  four  types  derive  from  three 
different  intervals  of  Ej  in  which  we  may  find  the  Jacobi  constant  of  the  orbit.  Table 
5-1  summarizes  these  four  types,  their  behaviors,  and  the  intervals  of  £7  in  which  they 
are  found. 

The  first  interval  only  has  an  upper  bound  to  Ej.  This  upper  bound  is  the  value 
of  Ej  at  the  Lagrange  points  Lj  and  L2.  Outside  the  bar  an  orbit  can  exist  with  an 
arbitrarily  low  value  of  £7  provided  that  it  is  located  at  a distance  sufficient  to  allow  it 
to  exist  above  the  zero  velocity  surface.  Within  the  bar,  however,  an  orbit  may  exist 
only  if  it  has  a Jacobi  constant  in  excess  of  the  minimum  £7  in  the  bar.  Also,  since  an 
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orbit  in  this  interval  of  Ej  has  insufficient  energy  to  traverse  the  L\  and  L2  points,  it 
remains  trapped  forever  either  inside  or  outside  the  bar  region.  Therefore,  this  interval 
of  values  of  Ej  leads  to  two  distinct  types  of  stochastic  orbits.  Figures  5-11  and  5-12 
give  examples  of  these  two  types  of  stochastic  behavior. 


Table  5-1:  Summary  of  basic  stochastic  orbit  types  and  their  behaviors. 


Interval  oi  Ej 

Type 

Distinct  behaviors 

(I)  Ej  < £y(Li) 

(1)  confined  to  bar 

circulation  in  bar 

(2)  confined  to  outer  disk 

circulation  in  outer  disk 

(0)  Ej{L{)  <Ej<  Ej{L^) 

(3)  confined  from  the 
vicinities  of  L4  and  L5,  but 
can  traverse  both  bar  and 
outer  disk 

(i)  circulation  in  bar 

(ii)  circulation  in  outer 
disk 

(iii)  circulation  around  L4 
or  L5 

(HI)  Ej  > Ej{L^) 

(4)  unconfined 
energetically 

circulation  everywhere 

kpc 


Figure  5-11:  Stochastic  orbit  trapped  within  the  bar.  The  value  of  its  Jacobi  constant 
{Ej  = -191307  km^  s"^)  is  slightly  less  than  that  of  the  Lagrange  points  L\  and 
L2  {Ej  = -191188  km^  s"^).  The  darker  circle  represents  corotation  at  5.5  kpc. 


138 


Figure  5-12:  Stochastic  orbit  confined  to  the  outer  disk.  The  value  of 
its  Jacobi  constant  (£y  = -194202  km^  s“^)  is  also  less  than  that  of  L\ 
and  La.  The  bar  and  spiral  potential  minima  are  drawn  for  reference. 

The  second  interval  of  £/  is  bounded  by  the  value  at  L\  (or  La)  and  the  value 
at  L4  (or  L5).  While  orbits  in  this  energy  range  cannot  move  freely  everywhere  (i.e. 
they  are  restricted  from  the  vicinities  of  L4  and  L5),  they  are  free  to  move  through 
the  bar  and  into  the  outer  disk  (Figure  5-13).  In  our  work  with  these  orbits  we  have 
noticed  three  subtypes  of  general  motion  (relative  to  the  rotating  frame)  which  they 
can  possess:  (1)  direct  circulation  within  the  bar,  (2)  circulation  around  the  Lagrange 
points  L4  and  L5  {direct  when  inside  these  points,  retrograde  when  outside),  and  (3) 
retrograde  circulation  in  the  outer  disk.  The  majority  of  these  orbits  also  form  loops 
which  involve  (temporary)  motion  counter  to  the  general  motion  (see  Figure  5-10).  We 
find  that,  given  sufficient  time,  all  of  the  orbits  in  this  range  of  £7  display  all  three 
subtypes  of  behavior;  however,  over  the  length  of  time  we  have  calculated  these  orbits 
in  our  models  (-1  billion  years),  they  may  exhibit  only  one  or  two  of  these  behaviors. 
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Figure  5-13:  Stochastic  orbit  which  traverses  both  the  bar  and  the  outer  disk.  The 
value  of  its  Jacobi  constant  {Ej  = -189692  km^  s“^)  is  greater  than  that 
of  L\  and  Li,  but  less  than  that  of  L4  and  L5  {Ej  = -187300  km^  s~^). 

These  stochastic  orbits  that  are  able  to  traverse  both  the  bar  and  disk  constitute  the 

so-called  “hot”  populations  observed  in  the  respective  2-D  and  3-D  N-body  simulations 

of  Sparke  and  Sellwood  (1987)  and  Pfenniger  and  Friedli  (1991).  As  can  be  seen 

in  Figure  5-13,  they  tend  to  make  many  loops  near  the  spiral  potential  minima  just 

beyond  the  ends  of  the  bar  (near  corotation),  and  we  find  that  this  “hot”  population 

is  instrumental  in  supporting  the  imposed  spiral  structure  in  this  region.  When  in 

the  region  of  the  bar,  however,  these  orbits  tend  to  remain  near  their  curves  of  zero 

velocity,  which  are  rather  broad  and  elliptical.  Hence,  the  same  orbits  that  support  spiral 

structure  may,  at  the  same  time,  weaken  strong,  thin  bar  structure.  For  comparison. 

Figure  5-14  shows  an  orbit  of  approximately  the  same  value  of  Ej  but  calculated  in  a 

spiral  perturbation  whose  amplitude  is  ten  times  greater.  While  we  see  that  the  orbit  still 

% 

traverses  the  bar  and  spiral,  it  spends  most  of  its  time  in  the  outer  disk  at  large  distances 


140 


and  does  not  support  the  imposed  bar  or  spiral  structure.  This  behavior,  therefore,  places 
upper  bounds  on  the  amplitudes  of  structures  which  might  be  enhanced  by  these  orbits. 


Figure  5-14:  Stochastic  orbit,  whose  Jacobi  constant  (Ej  = -189476  km  s“*) 
is  approximately  the  same  as  the  orbit  of  Figure  5-13,  calculated  in  a 
spiral  potential  of  ten  times  greater  amplitude  than  that  of  Figure  5-13. 

Finally,  we  have  the  type  of  stochastic  orbit  whose  Jacobi  constant  is  greater  than 
that  of  L4  (or  L5).  Figure  5-15  shows  an  orbit  of  this  type.  These  orbits  are  not 
restricted  from  any  region  energetically,  and  they  tend  to  fill  large  (almost  circular) 
regions  more  or  less  uniformly.  For  this  reason,  and  since  the  energies  of  these  orbits 
exceed  those  of  all  the  circular  orbits  in  the  axisymmetrized  disk,  we  assume  that  these 
orbits  are  negligibly  populated  and  ignore  them  in  our  models  as  being  of  only  very 
minor  importance.  In  the  following  section  we  give  estimates  of  the  numbers  of  each 

I 

type  of  orbit,  both  ordered  and  stochastic,  in  our  models. 
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Rgurc  5-15:  Stochastic  orbit  whose  Jacobi  constant  {Ej  - -183000 
km^  s"^)  exceeds  that  of  L4  and  L5  and  is  energetically  unconstrained. 

Proportions  of  Ordered  and  Stochastic  Orbits  in  the  Models ' 

Determining  the  relative  proportions  of  ordered  and  stochastic  orbits  in  our  models 
is  somewhat  difficult.  It  is  made  so  not  only  by  the  fact  that  we  do  not  keep  track  of 
whether  a given  dispersed  orbit  is  stochastic  or  ordered  in  the  modeling  process,  but 
also  by  the  fact  that  our  calculated  orbits  form  a discrete  representation  of  a continuous 
distribution  in  velocity  space  (i.e.  a particular  choice  of  a representative  dispersed  orbit 
may  lead  to  either  an  underemphasis  or  overemphasis  of  the  relative  importance  of 
stochasticity).  Other  techniques  of  self-consistent  modeling  do  not  suffer  from  this 
uncertainty.  For  example,  the  linear  programming  method  used  by  Schwarzschild 
(1979)  and  the  nonnegative  least-squares  method  employed  by  Pfenniger  (1984b)  require 
libraries  of  individual  orbits  from  which  their  self-consistent  models  are  constructed. 
Straightforward  combination  of  the  knowledge  of  the  regularity  or  stochasticity  of  each 
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orbit  with  its  population  weight  gives  directly  the  numbers  of  regular  versus  stochastic 
orbits.  Despite  these  difficulties,  we  have  devised  a method  by  which  we  can  roughly 
estimate  the  proportions  of  stochastic  versus  regular  orbits  in  our  models.  Moreover, 
we  can  estimate  how  the  stochastic  orbit  distribution  breaks  down  into  those  trapped  in 
the  bar,  those  confined  to  the  outer  disk,  and  those  in  the  “hot”  population,  respectively. 

The  data  at  our  disposal  for  this  estimation  process  are  of  two  types:  (1)  the  initial 
conditions  and  velocity  dispersions  of  the  model  orbit  families,  and  (2)  the  surfaces 
of  section  presented  earlier  in  this  chapter.  The  initial  conditions  of  an  orbit  family 
are  specified  by  radius  r and  radial  velocity  r along  the  x (intermediate  bar)  axis  as  a 
function  of  Jacobi  constant  (parameterized  by  the  circular  radius,  of  the  orbit  of  the 
same  Jacobi  constant  in  the  axisymmetric  case).  We  first  plot  r versus  Vc  for  each  family. 
The  region  of  trapping  about  a stable  periodic  family  corresponds  to  a range  of  values  of 
r.  We  estimate,  for  a given  value  of  r^,  the  amount  of  matter  trapped  on  regular  orbits 
around  this  family  by  integrating  our  assumed  Gaussian  velocity  distribution,  whose 
amplitude  and  width  are  given  by  the  weight  (i.e.  the  azimuthally  averaged  disk 
surface  density)  and  velocity  dispersion  a{r)  = cro  + rcOr  respectively,  over  this  trapping 
range.  Also,  if  there  are  more  than  one  family  at  a given  value  of  rc,  then  the  amplitude 
of  the  Gaussian  centered  on  the  wth  family  is  decreased  by  the  corresponding  relative 
weight  Wm  (see  Chapter  2).  Outside  the  trapping  ranges  and  up  to  the  limiting  curves 
the  contributions  of  the  Gaussians  are  assumed  to  be  stochastic.  Furthermore,  guided  by 
the  results  of  the  surfaces  of  section  (given  earlier  in  this  chapter),  we  assume  all  orbits 
dispersed  around  a “circular”  orbit  are  stochastic.  While  this  is  not  strictly  true,  the 
regular  orbits  that  do  exist  at  these  values  of  rc  have  very  large  radial  velocities  and  are 
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not  significantly  populated.  Finally,  we  break  down  the  stochastic  orbit  distributions 
further  by  sorting  them  according  to  the  value  of  their  Jacobi  constant.  Those  with 
Jacobi  constants  less  than  that  of  L\  and  Vc  inside  corotation  are  assumed  to  be  trapped 
in  the  bar.  Those  of  similar  Jacobi  constants  with  outside  corotation  are  assumed  to 
be  confined  to  the  disk.  Those  with  Jacobi  constants  greater  than  that  at  L\  are  taken 
to  be  the  “hot”  population. 

Table  5-2:  Estimated  breakdown  of  the  mass- weighted  orbit  populations  comprising  the  models 
of  Chapter  3 according  to  type  (trapped  versus  stochastic)  and  location  (bar,  disk,  or  both).  The 
errors  in  the  cited  figures  are  somewhat  uncertain,  but  are  estimated  to  be  of  the  order  of  5%. 


NGC  3992 


Bar  Population  (%) 

"Hot"  Population 
(%) 

Disk  Population 
(%) 

trapped 

74 

6 

stochastic 

5 

12 

3 

NGC  1073 

Bar  Population  (%) 

"Hot"  Population 

Disk  Population 

(%) 

(%) 

trapped 

63 

2 

stochastic 

14 

14 

7 

NGC  1398 

Bar  Population  (%) 

"Hot"  Population 

Disk  Population 

(%) 

(%) 

trapped 

59 

1 

stochastic 

10 

5 

25 

Table  5-2  gives  the  results  of  these  analyses  for  the  three  models  presented  in 
Chapter  3.  Several  conclusions  can  be  drawn  from  these  results.  First,  we  see  that  the 
bar  is  dominated  by  trapped  orbits  (primarily  around  the  x\  family).  For  example,  of  the 
total  number  of  orbits  confined  to  the  bar  region  in  the  best  model  of  NGC  3992,  fully 
94%  of  them  are  quasi-periodic.  This  agrees  with  the  results  of  several  other  studies,  for 
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example  Pfenniger  (1984b),  Teuben  and  Sanders  (1985),  Sparke  and  Sellwood  (1987), 
and  Pfenniger  and  Friedli  (1991).  A second  conclusion  is  that,  of  the  orbits  that  make 
up  the  outer  disk,  the  “hot”  population  comprises  a significant  fraction.  In  fact,  in  two 
of  the  models  (NGC  3992  and  NGC  1073)  this  group  is  larger  than  the  total  population 
confined  to  the  outer  disk  (regular  and  stochastic  orbits  combined).  In  their  study  of  a 
3-D  A-body  bar  and  disk  model,  Pfenniger  and  Friedli  (1991)  found  that  fully  31.3% 
of  the  particles  were  of  the  “hot”  variety  at  the  end  of  their  simulation  (55.9%  were 
located  in  the  bar  and  12.8%  were  in  the  disk).  Finally  we  note  that  there  seems  to 
be  slightly  too  many  estimated  orbits  in  the  bar  region.  For  example,  the  mass  of  the 
model  (bar  + exponential  disk)  inside  the  bar  radius  (5.5  kpc)  for  NGC  3992  is  62% 
of  the  total  mass  inside  r = 12.0  kpc.  From  Table  5-2  we  see  that  the  orbit  population 
estimates  give  some  79%.  Although  the  exact  cause  of  this  discrepancy  is  not  known, 
we  believe  that  it  is  probably  due  to  an  underestimation  of  the  “hot”  population.  It  is 
also  possible  that  the  discrepancy  is  due  simply  to  the  uncertainties  of  the  figures  in 
Table  5-2.  In  any  case  we  deem  the  two  previous  conclusions  to  be  firm. 


CHAPTER  6 
GAS  RESPONSE 


In  this  chapter  we  examine  the  response  of  gas  to  the  “best”  models  of  NGC  3992, 
NGC  1073,  and  NGC  1398.  In  order  to  calculate  the  gas  behavior  we  have  used  a 
two-dimensional  smoothed  particle  hydrodynamics  (SPH)  code.  Since  the  formalism 
of  our  SPH  approach  has  already  been  described  in  detail  in  Chapter  2,  we  restrict  our 
attention  here  to  the  details  of  the  simulations  and  to  the  presentation  and  interpretation 
of  the  results. 

There  are  two  basic  questions  to  be  answered  by  our  simulations:  (1)  does  a 
gaseous  component  exhibit  long-lived  spiral  structure  in  the  presence  of  the  strong 
stellar  spirals  of  our  “best”  models?  and  (2)  can  our  model  bars  alone  drive  spiral 
structure  in  the  gas?  The  answer  to  the  first  question  has  direct  bearing  on  the  overall 
success  of  our  models.  As  mentioned  in  Chapter  1,  the  spiral  arms  observed  in  most 
barred  galaxies  seem  to  exist  both  in  the  gaseous  component  and  in  the  underlying 
stellar  disk  (Elmegreen  and  Elmegreen  1985).  Thus  we  cannot  consider  as  ultimately 
successful  a model  in  which  the  gas  does  not  exhibit  long-lived  spiral  structure  along 
with  the  stars.  This  may  occur,  for  example,  if  the  stellar  spirals  are  so  strong  that  they 
either  disrupt  the  gaseous  features  soon  after  they  form  or  do  not  allow  them  to  form  at 
all.  The  second  question  addresses  the  issue  of  whether  stellar  spirals  are  necessary  at 
all  in  order  to  form  spiral  structure  in  the  gas.  Elmegreen  and  Elmegreen  (1985)  find 
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that  some  barred  galaxies  seem  to  possess  primarily  gaseous  arms.  This  leaves  open  the 
possibility,  mentioned  in  Chapter  1,  that  spiral  structure  is  solely  a bar-driven  feature 
of  the  gaseous  component.  With  these  two  questions  in  mind,  we  have  performed  the 
following  simulations. 

We  have  run  two  simulations  for  each  “best”  model.  The  first  employs  the  full 
amplitude  of  the  bar/spiral  perturbation  (in  addition  to  the  Ferrers  bar)  while  the  second 
employs  no  bar/spiral  perturbation  (only  the  Ferrers  bar).  In  all  simulations  the  initial 
gaseous  disk  is  represented  by  approximately  12,000  particles  placed  down  uniformly 
on  a Cartesian  grid.  The  radius  of  this  initial  particle  distribution  is  taken  to  be  roughly 
the  radius  of  the  outer  Lindblad  resonance  of  the  model.  The  particles  are  given  masses 
such  that  the  mass  of  the  whole  gaseous  disk  is  10%  of  the  mass  of  the  model  (out  to  the 
maximum  radius  of  the  initial  particle  distribution).  Also,  the  gas  particles  respond  only 
to  the  gravitational  forces  of  the  model  and  to  the  viscous  and  pressure  forces  among 
themselves  (i.e.  the  particles  are  not  self-gravitating).  The  initial  state  of  the  model  disk 
is  axisymmetric.  Flence,  the  amplitude  of  the  bar/spiral  perturbation  and  the  eccentricity 
of  the  Ferrers  bar  are  initially  set  to  zero.  The  particles  are  started  on  purely  circular 
trajectories  with  velocities  whose  magnitudes  are  determined  by  the  radial  forces  of  the 
initial  axisymmetric  model.  Over  the  first  rotation  of  the  pattern,  however,  the  bar/spiral 
amplitude  (if  used)  and  the  eccentricity  of  the  Ferrers  bar  are  allowed  to  increase  to 
their  full  values,  which  they  retain  for  the  duration  of  the  simulation. 

The  18  pages  of  figures  in  this  chapter  show  the  various  simulations  after  different 
numbers  of  rotations  of  the  bar  or  bar/spiral  pattern.  We  start  with  the  simulations  of 
NGC  3992,  and  proceed  to  those  of  NGC  1073  and  NGC  1398  in  turn.  Six  snapshots 
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per  simulation  are  given  showing  the  evolution  of  the  gas  structure.  For  a given  galaxy, 
the  first  six  snapshots  refer  to  the  run  with  the  bar/spiral  perturbation  imposed.  The 
second  six  snapshots  refer  to  the  run  with  the  Ferrers  bar  alone.  The  number  of  rotations 
of  the  bar  or  bar/spiral  pattern  is  noted  in  the  upper  right  hand  corner  of  the  figure.  In 
all  figures  the  frame  of  reference  is  the  one  corotating  with  the  bar/spiral  pattern  such 
that  the  bar  always  lies  along  the  y-axis  (i.e.  vertically). 

Figures  6-1  through  6-3  show  the  gas  response  in  the  model  of  NGC  3992  with 
the  bar/spiral  potential  imposed.  Well-defined  global  spiral  structure  is  seen  in  the  gas 
even  in  the  earliest  snapshot,  at  1.5  rotations  of  the  pattern  (top  panel  of  Figure  6-1).  At 
this  stage  the  gaseous  arms  exhibit  thin  ridges  of  high  density,  implying  the  existence 
of  shocks.  These  ridges  are  particularly  evident  in  the  outer  portions  of  the  disk,  from 
roughly  7 kpc  all  the  way  out  to  the  edge  of  the  disk  at  12  kpc.  We  also  notice  local 
enhancements  in  the  gas  density  near  the  ends  of  the  bar,  and  a slight  extension  of  the 
inner  termination  of  the  spirals  beyond  the  direction  of  the  bar,  giving  the  impression 
of  an  angular  offset  between  the  bar  and  spirals.  At  2.4  pattern  rotations  (bottom  panel 
of  Figure  6-1)  little  has  changed  except  for  the  beginnings  of  an  elongated  inner  ring 
surrounding  the  bar  and  a slight  tendency  for  the  gas  to  evacuate  the  regions  around  the 
Lagrange  points  L4  and  L5.  In  the  top  panel  (3.3  pattern  rotations)  and  bottom  panel 
(5.6  pattern  rotations)  of  Figure  6-2  we  see  the  continued  slow  development  of  the 
elongated  inner  ring  and  the  evacuation  of  the  regions  surrounding  L4  and  L5.  We  also 
notice  that  the  inner  spiral  structure,  extending  from  slightly  beyond  the  bar  radius  (-  6 
kpc)  to  about  9 kpc,  is  gradually  weakening.  In  fact,  the  outer  spiral  structure  (at  5.6 
rotations)  is  beginning  to  wrap  back  onto  itself  and  detach  from  the  weakening  inner 
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spiral  structure.  This  phenomenon  is  more  clearly  visible  in  the  panels  of  Figure  6-3, 
where  we  see  the  outer  spirals  forming  a rather  clumpy  ring  near  the  OLR  (10.6  kpc). 
Hiotelis  et  al.  (1993)  found  that  a similar  clumpy  gas  ring  formed  near  the  OLR  in 
their  SPH  calculation  of  the  response  of  viscous  gas  to  a purely  spiral-shaped  imposed 
potential.  Gaseous  spirals  still  exist  very  near  the  ends  of  the  bar,  but  in  between  the 
structure  is  significantly  weakened.  The  gas  structure  in  the  bar  is  essentially  unchanged 
from  the  earlier  snapshots.  We  note  that,  in  agreement  with  the  observations  of  NGC 
3992,  no  strong  offset  shocks  along  the  bar  are  present  at  any  time  during  the  calculation. 

A much  different  result  is  seen  in  the  calculation  where  no  outer  spiral  potential 
was  imposed  (Figure  6—4  through  Figure  6-6).  Only  weak  gaseous  spirals  are  observed 
at  any  time  during  this  run.  For  example,  at  1 .5  pattern  rotations  (top  panel  of  Figure 
6-4)  we  see  only  that  a short  gas  bar  has  formed,  and  that  gas  has  begun  to  accumulate 
somewhat  near  the  ends  of  the  bar  and  to  be  evacuated  from  regions  just  outside  the 
bar  (particularly  along  the  directions  of  the  bar  minor  axis).  Only  the  slightest  hint 
of  spiral  structure  is  visible.  At  2.5  rotations  (bottom  panel  of  Figure  6-4)  the  spiral 
structure  is  somewhat  stronger;  nevertheless,  the  spirals  that  do  appear  are  very  tightly 
wound  and  do  not  emanate  from  the  ends  of  the  bar  at  all.  In  fact,  these  spirals  are 
not  discernable  inside  a radius  of  approximately  7 kpc.  Figures  6-5  and  6-6  show  very 
similar  morphologies.  The  only  noticeable  change  in  the  gas  structure  from  3.5  pattern 
rotations  (top  panel  of  Figure  6-5)  to  10.6  pattern  rotations  (bottom  panel  of  Figure 
6-6)  is  the  slow  outward  movement  of  the  spirals  toward  the  edge  of  the  disk.  We 
do  notice  a tendency  for  the  gas  in  the  bar  to  accumulate  near  the  edges  of  the  bar 
by  10.6  pattern  rotations. 
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Figure  6- 
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: The  gas  response  in  the  best  model  of  NGC  3992  with  the  bar/spiral 
Top  panel:  1.5  pattern  rotations.  Bottom  panel:  2.4  pattern  rotations. 
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Figure  6-2:  The  gas  response  in  the  best  model  of  NGC  3992  with  the  bar/spiral  imposed 
(continued).  Top  panel:  3.3  pattern  rotations.  Bottom  panel:  5.6  pattern  rotations. 
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Figure  6-3:  The  gas  response  in  the  best  model  of  NGC  3992  with  the  bar/spiral  imposed 
(continued).  Top  panel:  8.1  pattern  rotations.  Bottom  panel:  10.4  pattern  rotations. 
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Figure  6-4:  The  gas  response  in  the  best  model  of  NGC  3992  with  only  the  Ferrers  bar 
present.  Top  panel:  1.5  pattern  rotations.  Bottom  panel:  2.5  pattern  rotations. 
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Figure  6-5:  The  gas  response  in  the  best  model  of  NGC  3992  with  only  the  Ferrers  bar 
present  (continued).  Top  panel:  3.5  pattern  rotations.  Bottom  panel:  5.3  pattern  rotations. 
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Figure  6-6:  The  gas  response  in  the  best  model  of  NGC  3992  with  only  the  Ferrers  bar 
present  (continued).  Top  panel:  7.6  pattern  rotations.  Bottom  panel:  10.6  pattern  rotations 
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Figure  6-7:  The  gas  response  in  the  best  model  of  NGC  1073  with  the  bar/spiral 
imposed.  Top  panel:  1.1  pattern  rotations.  Bottom  panel:  2.8  pattern  rotations. 
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Figure  6-8:  The  gas  response  in  the  best  model  of  NGC  1073  with  the  bar/spiral  imposed 
(continued).  Top  panel:  4.5  pattern  rotations.  Bottom  panel:  5.4  pattern  rotations. 
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Figure  6-9:  The  gas  response  in  the  best  model  of  NGC  1073  with  the  bar/spiral  imposed 
(continued).  Top  panel:  6.8  pattern  rotations.  Bottom  panel:  8.3  pattern  rotations. 
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Figure  6-10:  The  gas  response  in  the  best  model  of  NGC  1073  with  only  the  Ferrers  bar 
present.  Top  panel:  1.1  pattern  rotations.  Bottom  panel:  2.6  pattern  rotations. 
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Figure  6-11:  The  gas  response  in  the  best  model  of  NGC  1073  with  only  the  Ferrers  bar 
present  (continued).  Top  panel:  4.5  pattern  rotations.  Bottom  panel:  5.5  pattern  rotations. 
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Figure  6-12:  The  gas  response  in  the  best  model  of  NGC  1073  with  only  the  Ferrers  bar 
present  (continued).  Top  panel:  7.3  pattern  rotations.  Bottom  panel:  8.3  pattern  rotations. 
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Figure  6-13:  The  gas  response  in 
imposed.  Top  panel:  1.4  pattern 


the  best  model  of  NGC  1398  with  the  bar/spiral 
rotations.  Bottom  panel:  3.3  pattern  rotations. 
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Figure  6-14:  The  gas  response  in  the  best  model  of  NGC  1398  with  the  bar/spiral  imposed 
(continued).  Top  panel:  4.1  pattern  rotations.  Bottom  panel:  5.0  pattern  rotations. 
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Figure  6-15:  The  gas  response  in  the  best  model  of  NGC  1398  with  the  bar/spiral  imposed 
(continued).  Top  panel:  7.2  pattern  rotations.  Bottom  panel:  11.1  pattern  rotations. 
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Figure  6-16; 
present. 
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The  gas  response  in  the  best  model  of  NGC  1 398  with  only  the  Ferrers  bar 
Top  panel:  1.9  pattern  rotations.  Bottom  panel:  3.1  pattern  rotations. 
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Figure  6-17:  The  gas  response  in  the  best  model  of  NGC  1398  with  only  the  Ferrers  bar 
present  (continued).  Top  panel:  4.2  pattern  rotations.  Bottom  panel:  5.3  pattern  rotations. 
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Figure  6-18:  The  gas  response  in  the  best  model  of  NGC  1398  with  only  the  Ferrers  bar 
present  (continued).  Top  panel:  7.5  pattern  rotations.  Bottom  panel:  11.1  pattern  rotations 
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Figures  6-7  through  6-9  show  the  response  of  gas  to  the  model  of  NGC  1073 
with  the  imposed  bar/spiral  potential.  As  in  the  case  of  NGC  3992,  the  gas  displays 
well-defined  spiral  structure  from  the  earliest  snapshot,  at  1.1  rotations  of  the  pattern 
(top  panel  of  Figure  6-7).  At  this  early  stage  the  gas  density  is  enhanced  all  along  the 
spiral  arms  from  their  inward  terminations  near  the  ends  of  the  bar  to  the  edge  of  the 
gas  disk.  There  is,  however,  no  analogous  enhancement  of  the  bar  by  the  gas.  While  a 
small  amount  of  gas  has  congregated  at  the  center  of  the  galaxy  to  form  a very  short 
nuclear  bar  structure  there,  most  of  the  gas  originally  located  near  the  imposed  bar 
potential  has  already  begun  to  evacuate  the  region  by  1.1  pattern  rotations.  This  trend 
continues  until  a steady-state,  which  can  be  seen  in  the  next  snapshot  at  2.8  pattern 
rotations  (bottom  panel  of  Figure  6-7),  is  reached.  By  this  time  there  is  little  gas  in 
the  bar  and  essentially  none  near  the  Lagrange  points  L4  and  L5.  The  gas  that  does 
remain  near  the  bar  region  is  located  either  in  the  short  nuclear  bar  or  in  an  elongated 
inner  ring  surrounding  the  bar  region.  The  morphology  of  the  gas  changes  very  little 
and  very  slowly  from  this  point  (2.8  pattern  rotations)  until  the  end  of  the  calculations 
at  8.3  rotations  of  the  pattern  (bottom  panel  of  Figure  6-9).  There  is  a slight  tendency 
for  the  spiral  enhancements  to  weaken  near  the  ends  of  the  bar  and  to  become  smeared 
out  in  the  outer  disk.  These  two  effects  combine  to  make  the  gas  distribution  more 
and  more  like  an  outer  ring  with  time.  Here  also,  we  fail  to  note  the  development  of 
offset  shocks  along  the  bar  at  any  time. 

Figures  6-10  through  6-12  show  the  response  of  gas  to  the  model  of  NGC  1073 
without  the  imposed  bar/spiral  potential.  As  in  the  case  of  NGC  3992,  we  see  very 
different  gas  behavior  without  the  imposed  bar/spiral  potential.  At  1.1  pattern  rotations 
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(top  panel  of  Figure  6-10)  no  discernible  spiral  structure  is  present.  The  bar  region 
has  begun  to  be  evacuated  almost  completely  of  gas.  However,  as  in  the  case  of  the 
imposed  bar/spiral,  a small  quantity  of  gas  remains  at  the  center  and  forms  a nuclear 
bar.  By  2.6  pattern  rotations  (bottom  panel  of  Figure  6-10),  though,  we  do  see  some 
spiral  structure;  however,  it  is  quite  weak,  has  a very  small  pitch  angle,  and  does  not 
emanate  from  the  ends  of  the  bar.  No  spiral  structure  is  visible  inside  a radius  of  about 
4 kpc.  Figures  6-11  and  6-12  show  essentially  this  same  picture.  In  this  case  we  see 
no  tendency  for  the  gas  to  evolve  towards  an  outer  ring. 

Figures  6-13  to  6-15  depict  the  response  of  gas  to  the  model  of  NGC  1398  with 
the  bar/spiral  potential  imposed.  At  1.4  pattern  rotations  (top  panel  of  Figure  6—13) 
the  gas  exhibits  strong,  well-defined  spiral  structure.  As  in  the  analogous  run  for  the 
case  of  NGC  3992,  we  see  local  density  enhancements  in  the  gas  near  the  ends  of 
the  bar.  However,  by  3.3  pattern  rotations  (bottom  panel  of  Figure  6-13),  and  more 
clearly  by  4.1  rotations  (top  panel  of  Figure  6-14),  we  see  a dissolution  of  the  spiral 
structure  between  the  ends  of  the  bar  (-  4.8  kpc)  and  7 kpc.  At  the  same  time  the  gas 
that  forms  the  outer  spirals  (located  from  approximately  7 kpc  to  12  kpc)  is  disrupted, 
and  by  11.1  pattern  rotations  it  has  coalesced  into  a more  or  less  continuous  outer  ring 
located  roughly  at  the  OLR  (9.8  kpc).  Meanwhile,  the  gas  remaining  inside  corotation 
is  distributed  almost  uniformly  in  a thick  oval  bar.  Again,  we  do  not  observe  the 
formation  of  offset  shocks  along  the  bar. 

Figures  6-16  to  6-18  show  the  response  of  gas  to  this  same  model  of  NGC 
1398,  but  with  no  imposed  bar/spiral  potential.  In  this  case  we  see  an  evolution  very 
similar  to  the  case  of  NGC  3992  without  the  imposed  bar/spiral.  We  again  have  the 
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slight  enhancement  of  the  overall  gas  density  in  the  bar  at  early  times,  and  the  larger 
enhancements  near  the  edges  of  the  bar  at  later  times  (see  bottom  panel  of  Figure  6-18). 
Also,  there  is  a depletion  of  gas  in  the  regions  beyond  the  gas  bar,  especially  along  the 
directions  of  the  bar  minor  axis.  Outside  the  bar  we  have  the  generation  of  very  weak 
and  tightly  wound  spiral  waves  that  propagate  slowly  to  the  edge  of  the  disk  by  the  end 
of  the  calculation  (11.1  pattern  rotations,  bottom  panel  of  Figure  6-18). 

Examination  of  the  SPH  runs  yields  a couple  of  general  conclusions  concerning 
the  gas  response  in  our  models.  First,  in  cases  where  there  is  no  imposed  bar/spiral 
perturbation,  only  weak  spirals  are  found.  These  spirals  appear  rather  far  out  in  the 
outer  disk  and  soon  after  the  “formation”  of  the  Ferrers  bar  (i.e.  when  the  bar  reaches 
full  ellipticity),  propagate  slowly  toward  the  edge  of  the  gas  disk,  and  leave  the  rest  of 
the  disk  outside  the  bar  almost  completely  featureless  for  the  duration  of  the  run.  In  the 
cases  of  the  models  of  NGC  3992  and  NGC  1398,  a short,  stubby  gas  bar  forms  and 
the  remaining  area  within  the  bar  radius  becomes  somewhat  depleted  of  gas.  Thus  it 
seems  that  a bar  alone  is  insufficient  to  drive  strong  spirals  in  the  gas,  consistent  with 
the  findings  of  Ball  (1984,  1992),  Hunter  et  al.  (1988),  England  (1989)  and  England  et 
al.  (1990).  This  inability  to  drive  gaseous  spirals  is  due  to  the  fact  that  the  quadrupole 
field  of  the  Ferrers  bar  drops  off  too  rapidly  outside  the  bar. 

In  contrast,  though,  are  the  cases  where  the  bar/spiral  perturbation  is  imposed. 
Here  the  spiral  response  of  the  gas  is  very  strong.  In  the  case  of  NGC  1398,  in  fact,  the 
perturbation  is  so  strong  that  it  disrupts  the  gas  and  causes  it  to  form  an  outer  ring  after 
approximately  11  pattern  rotations  (Figure  6-15).  This  appears  to  be  consistent  with  the 
observed  HI  distribution  in  NGC  1398  (Moore,  private  communication).  Moreover,  the 
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gaseous  spiral  structure  just  outside  the  bar  is  substantially  weakened  and  has  just  about 
dissolved  long  before  this  (Figures  6-13  and  6-14).  The  optical  spiral  structure  in  NGC 
1398  is  indeed  quite  long,  thin,  and  continuous,  similar  to  the  panels  of  Figure  6-13. 
The  gas  response  in  the  case  of  NGC  3992,  though,  takes  the  form  of  relatively  long- 
lived  spirals  (Figures  6-1  through  6-3).  Therefore  it  appears  that  spiral  perturbations 
of  the  underlying  disk  are  necessary  in  order  to  generate  a long-lived  quasi-stationary 
response  in  the  gas. 

The  study  of  the  dynamics  of  gas  in  barred  galaxies  is  an  area  of  intense  ongoing 
research.  In  fact,  there  is  current  evidence  from  numerical  simulations  for  a dynamical 
coupling  of  the  stellar  and  gaseous  components  of  bars  (e.g.  Friedli  and  Benz  1992). 
However,  since  a full  consideration  of  this  aspect  of  our  models  would  take  us  too  far 
afield  from  the  present  topic  of  self-consistent  stellar  models,  we  defer  further  comment 
on  this  topic  to  future  investigations.  In  Chapter  7 we  summarize  the  conclusions  that 
we,  based  on  the  work  presented  in  this  thesis,  are  able  to  draw  concerning  the  nature 
of  self-consistent  stellar  models  of  barred  spiral  galaxies. 
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CHAPTER  7 
SUMMARY 


The  conclusions  of  this  dissertation  fall  into  three  basic  categories:  (1)  the  existence 
and  properties  of  self-consistent  stellar  models  of  barred  spiral  galaxies,  (2)  the  role  of 
stochastic  orbits  in  these  models,  and  (3)  the  response  of  gas  to  the  models.  We  end 
with  some  general  comments  regarding  our  models  and  directions  for  future  research. 

Self-Consistent  Models  of  Barred  Spirals 

Perhaps  the  most  important  conclusion  derived  from  this  work  is  that  self-consistent 
barred  spiral  models  exist.  Moreover,  they  exist  for  all  three  galaxies  considered  (i.e. 
for  a wide  range  of  Hubble  types).  The  parameters  of  the  most  successful  models 
reinforce  the  conclusions  of  previous  studies.  For  example,  we  find  in  each  case  that 
the  most  successful  model  places  corotation  near  the  end  of  the  bar.  Deviation  from  this 
condition  results  in  worse  self-consistency.  The  close  agreement  between  bar  length 
and  corotation,  as  mentioned  earlier,  is  also  found  in  most  A-body  simulations  and  orbit 
analyses.  Another  common  result  is  that  our  model  bars  are  made  up  almost  entirely  of 
regular,  elongated  orbits  trapped  around  the  Xi  family.  Also,  we  find  that  the  preferred 
phase  difference  between  the  bar  and  spirals  at  the  bar  radius  is  zero  (i.e.  the  arms 
connect  smoothly  to  the  ends  of  the  bar  with  no  angular  offset).  While  an  examination 
of  the  images  of  actual  barred  spiral  galaxies  seems  to  confirm  this  result,  Sellwood 
and  Sparke  (1988)  argue  that,  in  fact,  bars  and  spirals  represent  independent  patterns 
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with  different  angular  velocities.  The  obvious  objection  to  this  claim  is  that,  under 
this  assumption,  the  observed  phase  differences  between  the  ends  of  the  bars  and  the 
beginnings  of  the  spirals  should  be  randomly  distributed,  which  they  aren’t.  They  reply 
to  this  objection  by  noting  that  (1)  in  fact,  some  galaxies  do  exhibit  a phase  difference, 
and  (2)  the  bar  and  spirals  may  appear  to  the  eye  to  be  connected  when  actually  a 
phase  difference  exists.  The  numerical  simulations  of  Combes  and  Elmegreen  (1992), 
however,  support  the  idea  of  a single  pattern  speed.  In  only  one  case  did  a bar  and 
outer  spiral  pattern  exhibit  independent  angular  velocities.  Since  we  achieved  success 
using  uniform  pattern  speeds  for  the  bars  and  spirals,  we  were  not  forced  to  consider 
independent  pattern  speeds  and  therefore  did  not  check  this  possibility. 

While  a sample  of  three  is  too  small  to  allow  us  to  draw  firm  statistical  conclusions, 
several  properties  of  our  models  appear  to  correlate  with  Hubble  type.  First,  the 
magnitude  of  the  spiral  pitch  angle  is  larger  in  the  best  models  of  NGC  3992  and  NGC 
1073  than  in  the  best  model  of  NGC  1398.  In  addition,  and  somewhat  reassuringly, 
the  model  pitch  angles  agree  very  well  with  the  observed  pitch  angles.  Secondly,  the 
fact  that  the  Ferrers  bars  in  the  models  of  NGC  3992  and  NGC  1398  needed  to  be 
augmented  by  an  extension  of  the  spiral  perturbation  into  the  bar  region  indicates  that 
these  bars  are  qualitatively  different  than  the  bar  in  NGC  1073,  which  required  no 
such  augmentation.  This  qualitative  difference  seems  to  corroborate  the  findings  of 
Elmegreen  and  Elmegreen  (1985),  who  distinguish  between  flat  bars,  seen  typically  in 
early-type  galaxies,  and  exponential  bars,  predominantly  found  in  late-type  galaxies. 
While  these  authors  suggest  that  the  two  bar  types  result  from  different  formation 
mechanisms  and  have  different  extents  relative  to  their  corotation  radii,  we  find  for 
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these  two  different  types  successful  models  whose  bars  extend  all  the  way  to  corotation. 
Recent  work  (Combes  and  Elmegreen  1992)  has  given  support  to  the  idea  that  the 
exponential  bars  of  late-type  galaxies  may  not  be  limited  in  extent  by  corotation  but 
by  other  factors. 

The  Role  of  Stochastic  Orbits 

Another  important  result  of  our  work  is  that  stochastic  orbits  seem  to  play  a 
significant  role,  especially  near  and  just  beyond  the  ends  of  the  bar,  in  supporting  self- 
consistent  spiral  structure.  In  this  respect  the  most  important  stochastic  population  is 
the  so-called  “hot”  population,  consisting  of  orbits  with  sufficient  energy  to  be  able  to 
traverse  both  the  bar  and  the  outer  disk.  Particles  of  this  type  form  a significant  fraction 
of  the  total  in  the  respective  2-D  and  3-D  A-body  simulations  of  Sparke  and  Sellwood 
(1987)  and  Pfenniger  and  Friedli  (1991),  thus  we  expect  their  behavior  to  have  a large 
influence  on  long-lived  spiral  structure.  Stochastic  orbits  of  other  types  are  of  much 
less  important  in  our  models.  For  example,  stochastic  orbits  trapped  in  the  bar  region 
tend  to  fill  the  areas  delimited  by  their  zero  velocity  curves.  Since  these  areas  are  more 
oval  than  the  bar  itself,  stochastic  orbits  of  this  type  tend  to  weaken  the  bar.  These 
orbits,  however,  are  not  highly  populated  in  our  models  and  their  destructiveness  to  the 
bar  is  small.  Stochastic  orbits  confined  to  the  disk  tend  to  fill  rings  and  not  support 
imposed  spiral  structure.  Again  their  dynamical  role,  destructive  to  imposed  spirals,  is 
minimized  by  the  fact  that  they  are  rather  minimally  populated  both  relatively  (due  to 
increased  proportions  of  regular  orbits  as  the  spiral  amplitude  decreases  outward)  and 
absolutely  (due  to  decreased  matter  density  at  large  radii).  Finally,  stochastic  orbits 
with  energies  sufficient  to  allow  complete  exploration  of  the  bar  and  disk  regions  seem 
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also  to  be  unimportant  to  self-consistent  structure.  Not  only  do  they  not  support  the 
imposed  structure,  but  they  also  are  undoubtedly  less  populated  at  these  high  energy 
values.  A more  quantitative  exploration  of  the  latter  point  is  wanting,  however,  and 
could  (and  should)  be  included  in  the  quoted  results  of  future  A-body  simulations. 

The  Gas  Response 

A further  result,  derived  from  our  hydrodynamical  calculations,  is  that  bars  alone 
seem  incapable  of  driving  and  sustaining  gaseous  spiral  structure  in  the  outer  disks  of 
barred  spiral  galaxies.  The  presence  of  a source  of  nonaxisymmetric  forces  in  the  outer 
disk  appears  to  be  required.  Our  calculations  show  that  a spiral  stellar  wave  is  sufficient. 
This  result  agrees  well  with  the  observations  of  Elmegreen  and  Elmegreen  (1985),  who 
find  that  in  almost  all  cases  the  spiral  waves  of  barred  spirals  contain  enhancements 
both  in  the  young  objects  (0-B  stars,  HII  regions,  etc.,  which  are  tracers  of  the  gas) 
and  in  the  older  disk  stars.  Ball  (1984,  1992),  Hunter  et  al.  (1988),  England  (1989) 
and  England  et  al.  (1990)  have  found  that  oval  distortions  in  the  surface  density  of 
model  disks,  distortions  which  extend  beyond  the  bar  radius,  can  excite  and  maintain 
spiral  structure  in  the  gas  of  the  outer  disk.  However,  it  is  not  clear  at  this  time  whether 
these  oval  distortions  in  the  models  have  counterparts  in  real  galaxies.  Our  results  also 
indicate  that  self-consistent  stellar  spirals  may  be  so  strong  that  they  disrupt  the  gaseous 
spirals  (e.g.  in  the  gas  response  to  the  model  of  NGC  1398).  While  this  may  be  seen  as 
a weakness  of  our  models,  it  is  not  quite  clear  whether  these  results  actually  contradict 
the  observations.  For  example,  as  noted  in  Chapter  6,  the  gas  in  the  outer  disk  of  NGC 
1398  may  exist  primarily  in  an  outer  ring.  In  other  cases  (e.g.  NGC  3992)  the  gaseous 
spiral  structure  is  quite  long-lived. 
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Directions  for  Future  Research 

In  the  end  we  have  to  address  the  question  of  whether  our  models  represent  what 
is  actually  going  on  in  barred  spiral  galaxies.  Our  answer  is  a qualified  “yes.”  We 
have  obtained  self-consistent  models  (at  least  approximately)  of  our  program  galaxies, 
models  which  are  derived  from  and  consistent  with  the  observables  of  these  systems. 
The  properties  of  our  models  are  also  generally  consistent  with  the  results  of  both 
N-body  simulations  and  orbit  theory.  Still,  uncertainties  remain. 

How  stable  are  our  models?  This  question  has  not  been  addressed  at  all.  The 
assumption  of  quasi- stationary  bar  and  spiral  structure  is  what  has  allowed  us  to 
represent  the  bar  and  spiral  by  time-independent  potentials.  This  assumption  would 
be  invalid  if  the  evolution  of  the  system  is  rapid.  One  shortcoming  of  our  approach  is 
that  it  is  unable  to  answer  the  questions  of  stability  and  evolution,  and  we  must  resort 
to  other  means  (i.e.  fully  self-consistent  simulations)  for  these  answers. 

How  important  is  the  third  dimension?  Again,  this  question  has  not  been  fully 
addressed  in  our  models.  While  we  account  for  the  effect  that  varying  the  z scale 
height  of  the  bar/spiral  perturbation  has  on  the  resulting  response  surface  density,  the 
dynamics  which  generate  that  response  density  is  confined  to  the  two  dimensions  of  the 
principle  disk  plane.  This  simplification  does  not  introduce  any  serious  errors  as  long  as 
the  z thickness  is  small  compared  to  density  variations  in  the  plane.  If  this  is  the  case, 
the  z motion  decouples  from  the  motion  in  the  plane,  with  the  latter  determining  the 
structure  of  the  disk.  We  do,  however,  have  to  consider  the  z motions  in  the  case  of  bars, 
where  the  vertical  scale  height  is  comparable  to  the  scale  lengths  of  density  variations 
in  the  disk  plane.  In  this  case  the  z motion  may  couple  with  the  motion  in  the  plane 
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and  introduce  important  new  dynamical  phenomena  (e.g.  vertical  resonances,  complex 
instability,  Arnold  diffusion  of  stochastic  orbits,  etc.).  The  3-D  orbital  structure  of 
barred  galaxies  has  begun  to  be  investigated  in  recent  years,  most  notably  by  Pfenniger 
and  collaborators  (e.g.  Pfenniger  1984a,  1985;  Pfenniger  and  Friedli  1991),  but  by 
others  as  well  (e.g.  Patsis  and  Zachilas  1990;  Mahon  1992;  Zachilas  1993).  We 
need  to  generalize  our  approach  to  three  dimensions,  guided  by  the  results  of  these 
investigations. 

Finally,  how  does  gas  affect  the  dynamics  of  barred  galaxies?  We  have  presented 
in  this  thesis  the  passive  response  of  a gaseous  component  to  the  imposed  potentials 
of  our  models.  Surely  the  gas  must  play  a more  substantial  role  than  this.  There  is  a 
growing  body  of  evidence  (e.g.  Thomasson  et  al.  1990;  Moore  1992)  that  some  means 
of  holding  down  the  stellar  velocity  dispersion  is  required  in  order  to  obtain  long-lived 
stellar  spiral  structure.  Otherwise  the  stars  “heat  up”  and  spiral  structure  is  lost.  The 
dissipative  gas  component  may  yield  this  means.  We  know  that  the  gaseous  and  stellar 
components  are  coupled  via  star  formation  and  supernovae  explosions.  How  substantial 
is  this  coupling  and  what  are  its  effects  on  the  global  dynamics  of  the  system?  Also, 
what  is  the  role  of  gas  in  producing  the  nuclear  structures  (bars,  spirals)  sometimes 
seen  in  barred  galaxies?  All  of  these  questions  are  fundamental  to  our  understanding  of 
barred  spiral  galaxies,  and  it  is  toward  answering  these  questions  that  future  research 


should  be  directed. 


APPENDIX 

MODEL  BAR  QUANTITIES 


Here  is  summarized  the  quantities  associated  with  the  n = 2 Ferrers  bar  which  are 
needed  for  this  work. 

The  spatial  density  distribution  of  the  bar  is  given  by 

= (A-I) 

lO,  m>l, 
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= ^ + ^ + (a  > 6 > c > 0).  (A-2) 
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The  central  density  pc  is  related  to  the  total  mass  Mg  and  the  axes’  lengths  by 
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The  surface  density  projected  onto  the  z = 0 plane  is  gotten  by  integrating  the  spatial 
density  pix,y,z)  with  respect  to  z in  the  interval  where  the  spatial  density  is  nonzero. 
It  can  be  written  as 
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The  potential  in  the  z = 0 plane  is  given  by 
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axes  a,b,  and  c,  and  on  the  position  {x,y),  and  are  given  by 
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where  k)  and  k)  are  the  incomplete  elliptic  integrals,  sin  (j)  = 
^ = Y (a2Zc2)»  A is  the  unique  positive  solution  of 
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In  order  to  calculate  the  forces  due  to  the  bar,  the  partial  derivatives  of  the  potential 
with  respect  to  radius  and  azimuth  are  required: 
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Also,  the  second  partial  derivative  of  the  potential  with  respect  to  radius  is  given  by 
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